Open Case Studies - Rural and Urban Obesity

Motivation

BMI is used as a proxy for adiposity and is measured as an individual’s weight in kilograms divided by the individual’s height in meters squared. The different categories of weight are defined with the following BMI cutoffs as defined by the World Health Organization, the Centers for Disease Control and Prevention, and the National Institutes of Health:

BMI table source

The following chart can help you identify your BMI estimate:

BMI Chart source

Higher BMI (>30) is associated with higher rates of all-causes of mortality, as well as increased rates of type 2 diabetes, cancer, heart disease, and stroke.

An article published in Nature evaluated and compared the Body-Mass Index (BMI) of populations in rural and urban communities around the world:

NCD Risk Factor Collaboration (NCD-RisC). Rising rural body-mass index is the main driver of the global obesity epidemic in adults. Nature 569, 260–264 (2019).

This article challenged the widely-held view that increased urbanization was one of the major reasons for increased global obesity rates. This view came about because many countries around the world have shown increased urbanization levels in parallel with increased obesity rates. However this study demonstrated that this might not be the case; and that in fact for most regions around the world, BMI measurements are increasing in rural populations just as much, if not more so, than urban populations. This case study will evaluate the data reported in this article to explore regional and gender specific differences in the obesity rates around the world in 1985 and 2017. Most importantly we will test if there is a difference in obesity rates between rural and urban communities. To do this we will test if there is a difference between the means of the BMI measurements for each group.

Main Questions

Our main questions are:

1. Is there a difference between rural and urban BMI estimates around the world? How does this compare for males and females?
2. How have BMI estimates changed from 1985 to 2017?
3. How do different countries compare for BMI estimates- in particular, how does the United States compare to the rest of the world?

Learning Objectives

In this case study, we’ll walk you through importing data from a pdf, cleaning data, wrangling data, visualizing the data, and comparing two groups using well-established and commonly used packages, including stringr, tidyr, dplyr, purrr, and ggplot2. We will especially focus on using packages and functions from the Tidyverse. The Tidyverse is a library of packages created by the chief scientist at RStudio, Hadley Wickham. While some students may be familiar with previous R programming packages, these packages make data science in R especially efficient.

We will begin by loading the packages that we will need:

library(here)
library(pdftools)
library(stringr)
library(readr)
library(dplyr)
library(tibble)
library(magrittr)
library(glue)
library(purrr)
library(tidyr)
library(ggplot2)
library(ggrepel)
library(cowplot)

Package	Use
here	to easily load and save data
pdftools	to read a pdf into R
stringr	to manipulate the text within the pdf of the data
readr	to manipulate the text within the pdf of the data into individual lines
dplyr	to arrange/filter/select subsets of the data
tibble	to create data objects that we can manipulate with dplyr/stringr/tidyr/purrr
magrittr	to use the %<>% pipping operator
glue	to paste or combine character strings and data together
purrr	to perform functions on all columns of a tibble
tidyr	to convert data from wide to long format
ggplot2	to make visualizations with multiple layers
ggrepel	to allow labels in figures not to overlap
cowplot	to allow plots to be combined

---

The first time we use a function, we will use the :: to indicate which package we are using. Unless we have overlapping function names, this is not necessary, but we will include it here to be informative about where the functions we will use come from.

Context

The measurement of BMI has some limitations that are well recognized, as it does not account for the composition of body mass, the location of body fat, or the contribution of body frame size. However, BMI has been a useful health indicator for risk for many diseases and conditions particularly when combined with other risk factor information.

What are the data?

We will be using data located within a table of the supplementary material for the NCD-RisC paper referenced above. This is a pdf that can be found freely available online.

Here you can see that the data contains mean BMI values for both men and women in various countries at the national level, as well as the mean BMI values for the rural and urban areas of these countries for both 1985 and 2017.

The data within the parentheses are the 95 % credible interval (CIs) ranges for the mean BMI estimates. The authors provide these CIs as a guide to understand how likely the estimate is for the true population mean BMI. A wider range suggests that the estimate is less accurate, as there are more possible values for the true mean with credible evidence.

Note: While gender and sex are not actually binary, the data presented that is used in this analysis only contains data for groups of individuals described as men or women.

Data Import

First let’s download the data:

url <- "https://static-content.springer.com/esm/art%3A10.1038%2Fs41586-019-1171-x/MediaObjects/41586_2019_1171_MOESM1_ESM.pdf"
utils::download.file(url, "paper_supplement.pdf")

Now that we have the obesity pdf, we will read it in to R using the pdftools package:

pdf_obesity<-pdftools::pdf_text(here("paper_supplement.pdf"))
# PDF error: 
# "Expected the optional content group list, 
# but wasn't able to find it, or it isn't an Array"

# Errors like this can be common with PDFs
# We can take a look at the data to see if it still looks as expected

Let’s take a look at the data- the summary() function helps us to look at the structure of R objects.

summary(pdf_obesity)

   Length     Class      Mode 
       63 character character 

We can see that we have 63 elements that are character strings. You may also notice that the original PDF has 63 pages. Let’s take a look at some of these elements.

pdf_obesity[1] # this looks like the first page

[1] "Letter                                                                                https://doi.org/10.1038/s41586-019-1171-x\nRising rural body-mass index is the main driver of\nthe global obesity epidemic in adults\nNCD Risk Factor Collaboration (NCD-RisC)*\n*A list of authors and their affiliations appears in the online version of the paper.\n2 6 0 | N A T U RE | V O L 5 6 9 | 9 M A Y 2 0 1 9\n"

pdf_obesity[2] # this looks like the second page

[1] "Supplementary Information. Statistical model for estimating BMI trends by rural and urban\nplace of residence.\n                                          1\n"

pdf_obesity[55] 
# this looks like data in a the table we want

We can see that the output looks pretty similar to the pages of the pdf, but the spacing is a bit awkward. Note that the way the data is displayed is partially influenced by the width setting of the RStudio window.

We are interested in a supplementary Table 3. which has multiple pages and includes the same header on each page; we can use that to determine what elements of our pdf_obesity character strings include our table. We will use the str_detect() function of the stringr package to search for the elements that contain text that is consistently in the header. The output of of this function will show which elements of the object (in this case pages of the pdf) include this pattern indicated as a “TRUE” or “FALSE”.

stringr::str_detect(pattern = "Age-standardised mean BMI", string = pdf_obesity)

 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[23] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[34] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[45] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE
[56]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE

# "Age-standardised mean BMI" is part of the header in the table on every page.
# This shows the pages that contain our table of interest (as TRUE values).

Now we will extract just the data for the table now and call it rural_urban. To do this we will use the str_subset() function of the stringr package.

rural_urban <-stringr::str_subset(pattern ="Age-standardised mean BMI", 
                                  string =pdf_obesity)

str(rural_urban) 

 chr [1:10] "                                                                          2                                    "| __truncated__ ...

# The str() function is similar to the summary function
# There are 10 pages worth of elements in our rural_urban object

Let’s check the first and last page:

cat(rural_urban[1]) 
# Using the cat() function which stands for Concatenate And Print 
# This will allow for the "\n" values to be shown as spaces

This looks the same as the beginning… how about the end?

cat(rural_urban[10])

Great! Our rural_urban object looks like it contains the entire Supplementary 3 table, as both the beginning and the end include the data we expected.

Data Wrangling

At this point we have large strings now for each page of the table, but this is not very convenient to work with. Now we will wrangle the data into a more usable form. Ideally we would like our data to be in some sort of tabular form.

Separate the data into lines

First it would be useful to separate the data into lines or rows of the table.

rural_urban <-readr::read_lines(rural_urban)
str(rural_urban) 

 chr [1:461] "                                                                          2                                    "| __truncated__ ...

# now we have 461 lines
rural_urban[6] 

[1] "Afghanistan          Women   20.6 (18.4-22.8)      20.1 (17.8-22.4)     23.2 (20.9-25.4) 24.4 (23.3-25.4)      23.6 (22.5-24.8)     26.3 (25.1-27.4)"

# we can see that line 6 shows the data for females in Afghanistan  

Removing excess white-space

We also have a lot of white-space… let’s get rid of the excess white spaces using str_squish() also from the stringr package.

rural_urban<-stringr::str_squish(rural_urban) 
head(rural_urban)

[1] "2 2"                                                                                                                    
[2] "Age-standardised mean BMI in 1985 (kg/m ) Age-standardised mean BMI in 2017 (kg/m )"                                    
[3] "Country Sex"                                                                                                            
[4] "National Rural Urban National Rural Urban"                                                                              
[5] "Men 20.2 (17.8-22.7) 19.7 (17.2-22.2) 22.4 (20.0-25.0) 22.8 (20.3-25.3) 22.5 (20.0-25.0) 23.6 (21.0-26.1)"              
[6] "Afghanistan Women 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (22.5-24.8) 26.3 (25.1-27.4)"

Now it is much easier to see the data.

If we look at the end of the first page of the table and the start of the second we can see that the header information is repeated, as well as a line with the page number and an empty line, and a line that says “2 2”.

rural_urban[44:56]

 [1] "Benin Women 20.6 (19.1-22.0) 20.1 (18.5-21.6) 21.8 (20.3-23.3) 24.3 (23.5-25.0) 23.2 (22.5-24.0) 25.5 (24.7-26.3)"  
 [2] "Men 24.3 (21.8-26.8) na* 24.3 (21.8-26.8) 26.7 (24.3-29.2) na* 26.7 (24.3-29.2)"                                    
 [3] "Bermuda Women 25.4 (22.2-28.6) na* 25.4 (22.2-28.6) 28.5 (25.3-31.6) na* 28.5 (25.3-31.6)"                          
 [4] "Men 20.6 (19.1-22.1) 20.3 (18.7-21.8) 23.1 (21.5-24.6) 23.5 (22.8-24.3) 23.1 (22.2-23.9) 24.2 (23.3-25.2)"          
 [5] "Bhutan Women 20.7 (18.4-22.9) 20.3 (18.0-22.6) 23.0 (20.8-25.3) 24.6 (23.7-25.5) 23.8 (22.7-24.9) 25.9 (24.7-27.0)" 
 [6] "51"                                                                                                                 
 [7] ""                                                                                                                   
 [8] "2 2"                                                                                                                
 [9] "Age-standardised mean BMI in 1985 (kg/m ) Age-standardised mean BMI in 2017 (kg/m )"                                
[10] "Country Sex"                                                                                                        
[11] "National Rural Urban National Rural Urban"                                                                          
[12] "Men 23.3 (21.0-25.5) 22.9 (20.6-25.2) 23.6 (21.3-25.9) 26.1 (23.9-28.4) 25.3 (23.1-27.5) 26.5 (24.2-28.8)"          
[13] "Bolivia Women 23.8 (22.3-25.3) 23.0 (21.4-24.5) 24.6 (23.1-26.1) 27.9 (26.5-29.3) 27.0 (25.6-28.5) 28.3 (26.8-29.7)"

Although the header was necessary on all of the pages of the pdf version of the table, we only need that information once in our data.

Removing unnecessary repeated header information

So, let’s remove all the header information and the page number lines from the rural_urban object, then we will make a single line header for the beginning. One way to do this is to find all lines that include either “women” or “men” and only keep this data. First let’s see how many lines include “women” or “men”. We can do this by first identifying the lines that contain these patterns and then check the length of the resulting vector for the line numbers that contain these patterns.

stringr::str_which(string = rural_urban, pattern = "Women|Men")

  [1]   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21
 [18]  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38
 [35]  39  40  41  42  43  44  45  46  47  48  55  56  57  58  59  60  61
 [52]  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78
 [69]  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95
 [86]  96  97  98 105 106 107 108 109 110 111 112 113 114 115 116 117 118
[103] 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
[120] 136 137 138 139 140 141 142 143 144 145 146 147 148 155 156 157 158
[137] 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175
[154] 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192
[171] 193 194 195 196 197 198 205 206 207 208 209 210 211 212 213 214 215
[188] 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232
[205] 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 255
[222] 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272
[239] 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289
[256] 290 291 292 293 294 295 296 297 298 305 306 307 308 309 310 311 312
[273] 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329
[290] 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346
[307] 347 348 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369
[324] 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386
[341] 387 388 389 390 391 392 393 394 395 396 397 398 405 406 407 408 409
[358] 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426
[375] 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443
[392] 444 445 446 447 448 455 456 457 458

length(stringr::str_which(string = rural_urban, pattern = "Women|Men"))

[1] 400

# Here the "|" indicates that we are looking for either  "Women" or "Men"

OK, so this looks correct. This includes most lines but there are gaps where are header is located. It looks like there are 400 lines in our table that aren’t headers.

rural_urban <- str_subset(string = rural_urban, pattern = "Women|Men")

We can check our data now using either head() or glimpse().

dplyr::glimpse(rural_urban)

 chr [1:400] "Men 20.2 (17.8-22.7) 19.7 (17.2-22.2) 22.4 (20.0-25.0) 22.8 (20.3-25.3) 22.5 (20.0-25.0) 23.6 (21.0-26.1)" ...

Great, now our rural_urban object 400 lines instead of 461, like it did before when it contained the redundant header information.

head(rural_urban)

[1] "Men 20.2 (17.8-22.7) 19.7 (17.2-22.2) 22.4 (20.0-25.0) 22.8 (20.3-25.3) 22.5 (20.0-25.0) 23.6 (21.0-26.1)"              
[2] "Afghanistan Women 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (22.5-24.8) 26.3 (25.1-27.4)"
[3] "Men 25.2 (23.9-26.5) 25.0 (23.7-26.4) 25.4 (24.0-26.7) 27.0 (26.0-27.9) 26.9 (25.9-27.9) 27.0 (26.0-28.0)"              
[4] "Albania Women 26.0 (24.1-27.9) 26.1 (24.1-28.1) 25.9 (23.9-27.8) 26.0 (24.8-27.2) 26.2 (24.8-27.5) 25.9 (24.6-27.2)"    
[5] "Men 22.1 (20.8-23.3) 21.8 (20.5-23.1) 22.3 (21.0-23.6) 25.1 (24.5-25.7) 24.8 (24.1-25.4) 25.2 (24.6-25.9)"              
[6] "Algeria Women 24.0 (22.2-25.7) 23.3 (21.4-25.1) 24.8 (22.9-26.6) 27.4 (26.7-28.0) 27.0 (26.3-27.8) 27.5 (26.7-28.2)"    

The head() function shows us the first rows or lines of the data, while the glimpse() function provides us information about the total size of the object and shows us the first line or row.

Great! So now our data looks much better but we need to add back our header and we would like this to only be a single line to make it easy to transform our data into a table or table-like object.

Here is what our header used to look like:

Dealing with spacing

First let’s try splitting our header-less data into columns based on spaces using the str_split() function, where we specify that we are splitting the data based on the pattern of a space (the space is included in quotes):

(str_split(string = rural_urban, pattern = " ", simplify = TRUE))[1:10,] 

      [,1]          [,2]    [,3]          [,4]          [,5]         
 [1,] "Men"         "20.2"  "(17.8-22.7)" "19.7"        "(17.2-22.2)"
 [2,] "Afghanistan" "Women" "20.6"        "(18.4-22.8)" "20.1"       
 [3,] "Men"         "25.2"  "(23.9-26.5)" "25.0"        "(23.7-26.4)"
 [4,] "Albania"     "Women" "26.0"        "(24.1-27.9)" "26.1"       
 [5,] "Men"         "22.1"  "(20.8-23.3)" "21.8"        "(20.5-23.1)"
 [6,] "Algeria"     "Women" "24.0"        "(22.2-25.7)" "23.3"       
 [7,] "Men"         "33.7"  "(32.7-34.7)" "32.6"        "(31.7-33.5)"
 [8,] "American"    "Samoa" "Women"       "34.3"        "(33.1-35.6)"
 [9,] "Men"         "25.0"  "(22.5-27.4)" "25.3"        "(22.8-27.7)"
[10,] "Andorra"     "Women" "25.2"        "(22.0-28.4)" "25.4"       
      [,6]          [,7]          [,8]          [,9]         
 [1,] "22.4"        "(20.0-25.0)" "22.8"        "(20.3-25.3)"
 [2,] "(17.8-22.4)" "23.2"        "(20.9-25.4)" "24.4"       
 [3,] "25.4"        "(24.0-26.7)" "27.0"        "(26.0-27.9)"
 [4,] "(24.1-28.1)" "25.9"        "(23.9-27.8)" "26.0"       
 [5,] "22.3"        "(21.0-23.6)" "25.1"        "(24.5-25.7)"
 [6,] "(21.4-25.1)" "24.8"        "(22.9-26.6)" "27.4"       
 [7,] "34.0"        "(32.9-35.1)" "34.3"        "(33.0-35.6)"
 [8,] "34.1"        "(33.0-35.2)" "34.4"        "(33.0-35.8)"
 [9,] "25.0"        "(22.5-27.4)" "26.8"        "(24.4-29.2)"
[10,] "(22.2-28.7)" "25.2"        "(22.0-28.4)" "25.3"       
      [,10]         [,11]         [,12]         [,13]        
 [1,] "22.5"        "(20.0-25.0)" "23.6"        "(21.0-26.1)"
 [2,] "(23.3-25.4)" "23.6"        "(22.5-24.8)" "26.3"       
 [3,] "26.9"        "(25.9-27.9)" "27.0"        "(26.0-28.0)"
 [4,] "(24.8-27.2)" "26.2"        "(24.8-27.5)" "25.9"       
 [5,] "24.8"        "(24.1-25.4)" "25.2"        "(24.6-25.9)"
 [6,] "(26.7-28.0)" "27.0"        "(26.3-27.8)" "27.5"       
 [7,] "34.6"        "(33.1-35.9)" "34.2"        "(32.9-35.6)"
 [8,] "35.3"        "(33.7-36.9)" "35.0"        "(33.1-36.9)"
 [9,] "26.8"        "(24.3-29.2)" "26.8"        "(24.3-29.3)"
[10,] "(22.1-28.6)" "25.2"        "(21.9-28.5)" "25.3"       
      [,14]         [,15]         [,16] [,17] [,18]
 [1,] ""            ""            ""    ""    ""   
 [2,] "(25.1-27.4)" ""            ""    ""    ""   
 [3,] ""            ""            ""    ""    ""   
 [4,] "(24.6-27.2)" ""            ""    ""    ""   
 [5,] ""            ""            ""    ""    ""   
 [6,] "(26.7-28.2)" ""            ""    ""    ""   
 [7,] ""            ""            ""    ""    ""   
 [8,] "35.4"        "(33.7-37.1)" ""    ""    ""   
 [9,] ""            ""            ""    ""    ""   
[10,] "(22.1-28.6)" ""            ""    ""    ""   

# Here we will take a look at just the first 10 lines using the [1:10,]

This almost worked, but unfortunately country names that have spaces will be a problem. We can see that American Samoa has been divided into two columns and all subsequent columns are shifted.

So now let’s try to extract the country information, by separating the country information from the sex information when the sex is female.

Sex always starts with either a capital “W” if the gender is female. #We need to use a space before the “W” otherwise we will split some of the country names if the names starts with “W”.

Here we will also introduce the concept of piping, which uses the %>%. This is really useful when we have multiple steps, which we will show soon.

So first we will select just the data for Women and then split this data based on the pattern " Women" (Note that the space is included - so that it is not within the country data).

Women <- str_subset(string = rural_urban, pattern = "Women") 
# First let's subset the data to only the lines that contain "Women"

country_split <-Women %>%
  stringr::str_split(pattern= " Women") %>%
  unlist() 


#we can do the same thing with two pipes:
# first we select the "Women" data and then we split the data by " Women"
# finally we assign the output of both steps to the object country_split

country_split <- str_subset(string = rural_urban, pattern = "Women") %>%
  stringr::str_split(pattern= " Women") %>%
  unlist()

head(country_split)

[1] "Afghanistan"                                                                                           
[2] " 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (22.5-24.8) 26.3 (25.1-27.4)"
[3] "Albania"                                                                                               
[4] " 26.0 (24.1-27.9) 26.1 (24.1-28.1) 25.9 (23.9-27.8) 26.0 (24.8-27.2) 26.2 (24.8-27.5) 25.9 (24.6-27.2)"
[5] "Algeria"                                                                                               
[6] " 24.0 (22.2-25.7) 23.3 (21.4-25.1) 24.8 (22.9-26.6) 27.4 (26.7-28.0) 27.0 (26.3-27.8) 27.5 (26.7-28.2)"

Now we can see that Country is always the odd rows and the rest of the data is the rest of the rows. We can select the odd rows using the Modulus (Remainder from division) operator. This operator looks like this %%. All odd values such as 3, 7, or 15 when divided by 2 would leave a remainder of 1. We will select just these rows using the filter() function of dplyr.

country <-tibble(country_split) %>% 
  dplyr::filter(row_number() %% 2 == 1) 

We can take a look to make sure that all the country names look as expected:

country

# A tibble: 200 x 1
   country_split      
   <chr>              
 1 Afghanistan        
 2 Albania            
 3 Algeria            
 4 American Samoa     
 5 Andorra            
 6 Angola             
 7 Antigua and Barbuda
 8 Argentina          
 9 Armenia            
10 Australia          
# … with 190 more rows

Looks good!

The even rows are the rows with the data for women. These row values have a remainder of 0 when divided by 2.

Women_BMI <-tibble(country_split) %>% 
  dplyr::filter(row_number() %% 2 == 0)
head(Women_BMI)

# A tibble: 6 x 1
  country_split                                                            
  <chr>                                                                    
1 " 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23…
2 " 26.0 (24.1-27.9) 26.1 (24.1-28.1) 25.9 (23.9-27.8) 26.0 (24.8-27.2) 26…
3 " 24.0 (22.2-25.7) 23.3 (21.4-25.1) 24.8 (22.9-26.6) 27.4 (26.7-28.0) 27…
4 " 34.3 (33.1-35.6) 34.1 (33.0-35.2) 34.4 (33.0-35.8) 35.3 (33.7-36.9) 35…
5 " 25.2 (22.0-28.4) 25.4 (22.2-28.7) 25.2 (22.0-28.4) 25.3 (22.1-28.6) 25…
6 " 21.3 (18.0-24.6) 20.9 (17.6-24.3) 22.7 (19.3-26.0) 24.4 (21.2-27.7) 23…

Great! Now we have a list of the countries that can be used for both the male and female data.

It’s always a good idea to check that your data objects are the size you expect when wrangling. We can do so with the dim() function, which shows us the dimensions of data objects.

dim(Women_BMI)

[1] 200   1

Great! There are 200 rows like we expected.

We do however need to add " Women" back to this data.

To do this we will use the glue() function of the glue package.

Women <-glue::glue("Women{Women_BMI$country_split}")
head(Women)

Women 20.6 (18.4-22.8) 20.1 (17.8-22.4) 23.2 (20.9-25.4) 24.4 (23.3-25.4) 23.6 (22.5-24.8) 26.3 (25.1-27.4)
Women 26.0 (24.1-27.9) 26.1 (24.1-28.1) 25.9 (23.9-27.8) 26.0 (24.8-27.2) 26.2 (24.8-27.5) 25.9 (24.6-27.2)
Women 24.0 (22.2-25.7) 23.3 (21.4-25.1) 24.8 (22.9-26.6) 27.4 (26.7-28.0) 27.0 (26.3-27.8) 27.5 (26.7-28.2)
Women 34.3 (33.1-35.6) 34.1 (33.0-35.2) 34.4 (33.0-35.8) 35.3 (33.7-36.9) 35.0 (33.1-36.9) 35.4 (33.7-37.1)
Women 25.2 (22.0-28.4) 25.4 (22.2-28.7) 25.2 (22.0-28.4) 25.3 (22.1-28.6) 25.2 (21.9-28.5) 25.3 (22.1-28.6)
Women 21.3 (18.0-24.6) 20.9 (17.6-24.3) 22.7 (19.3-26.0) 24.4 (21.2-27.7) 23.3 (20.0-26.7) 25.7 (22.4-29.1)

# here we will paste " Women" in fromt of the Women_BMI data
# note that there will be a space in between
# This is because the space is already in front of the first BMI value

Let’s grab the male data:

# remember our rural_urban object contains male data for the odd rows

Men <-tibble(rural_urban) %>% 
  dplyr::filter(row_number() %% 2 == 1) 
head(Men$rural_urban)

[1] "Men 20.2 (17.8-22.7) 19.7 (17.2-22.2) 22.4 (20.0-25.0) 22.8 (20.3-25.3) 22.5 (20.0-25.0) 23.6 (21.0-26.1)"
[2] "Men 25.2 (23.9-26.5) 25.0 (23.7-26.4) 25.4 (24.0-26.7) 27.0 (26.0-27.9) 26.9 (25.9-27.9) 27.0 (26.0-28.0)"
[3] "Men 22.1 (20.8-23.3) 21.8 (20.5-23.1) 22.3 (21.0-23.6) 25.1 (24.5-25.7) 24.8 (24.1-25.4) 25.2 (24.6-25.9)"
[4] "Men 33.7 (32.7-34.7) 32.6 (31.7-33.5) 34.0 (32.9-35.1) 34.3 (33.0-35.6) 34.6 (33.1-35.9) 34.2 (32.9-35.6)"
[5] "Men 25.0 (22.5-27.4) 25.3 (22.8-27.7) 25.0 (22.5-27.4) 26.8 (24.4-29.2) 26.8 (24.3-29.2) 26.8 (24.3-29.3)"
[6] "Men 20.5 (17.9-23.1) 20.2 (17.6-22.9) 21.4 (18.8-24.0) 22.6 (20.0-25.1) 22.0 (19.4-24.6) 23.2 (20.6-25.9)"

# Again let's make sure we have the correct number of rows 

dim(Men) 

[1] 200   1

How about our number of columns?

If we try splitting our data by space again, will it have the expected number of columns? What about the rows that contain na* values?

# Let's just take the Men data that contains na* values
# This column is called rural_urban
unknown <-Men %>%filter(str_detect(pattern ="na\\*", string = rural_urban)) 

# Now we can try splitting by a space
tibble::as_tibble(str_split(unknown$rural_urban, " ", simplify = TRUE))

# A tibble: 5 x 12
  V1    V2    V3     V4    V5     V6    V7    V8    V9    V10   V11   V12  
  <chr> <chr> <chr>  <chr> <chr>  <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 Men   24.3  (21.8… na*   24.3   (21.… 26.7  (24.… na*   26.7  (24.… ""   
2 Men   22.4  (21.4… 21.5  (20.2… 22.4  (21.… 24.8  (23.… na*   24.8  (23.…
3 Men   32.6  (32.0… na*   32.6   (32.… 32.9  (31.… na*   32.9  (31.… ""   
4 Men   22.8  (22.2… na*   22.8   (22.… 24.4  (23.… na*   24.4  (23.… ""   
5 Men   29.6  (28.1… 29.6  (28.1… na*   32.3  (31.… 32.3  (31.… na*   ""   

Dealing with NA values

So close! Notice that the "na*" values have shifted the subsequent values within the columns because typically there is a space between the BMI and the credible intervals. Here we can see this data in our original pdf:

We need to replace our na* values with something that includes a space so that when we separate our data by space we will have two values instead of one when we have an na* . Therefore, na* na* should work.

class(Men$rural_urban)

[1] "character"

# In order to use the functions within the stringr package, 
# our Men data needs to be of character class
Men<-as.character(Men$rural_urban) 
# This function allows us to change the factor data 
# within Men$rural_urban to the character class
# We are also changing the structure so that the data 
# is not within a vector called rural_urban
Men<-stringr::str_replace_all(string = Men, 
                              pattern = "na\\*", 
                              replacement = "na\\* na\\*") 
# The \\ are necessary because the * is a special character 
# The * would typically indicate any possible value, 
# but here we actually want a "*" instead
# Thus the double backslash does that for us
# Here we are replacing all occurences of the na* values 
#(thus str_replace_all instead of str_replace) with na* na*. 

Let’s check that it worked…

Men[20:30]

 [1] "Men 20.7 (19.4-22.1) 20.3 (18.9-21.7) 21.7 (20.3-23.1) 22.7 (22.0-23.3) 22.1 (21.3-22.8) 23.4 (22.7-24.1)"
 [2] "Men 24.3 (21.8-26.8) na* na* 24.3 (21.8-26.8) 26.7 (24.3-29.2) na* na* 26.7 (24.3-29.2)"                  
 [3] "Men 20.6 (19.1-22.1) 20.3 (18.7-21.8) 23.1 (21.5-24.6) 23.5 (22.8-24.3) 23.1 (22.2-23.9) 24.2 (23.3-25.2)"
 [4] "Men 23.3 (21.0-25.5) 22.9 (20.6-25.2) 23.6 (21.3-25.9) 26.1 (23.9-28.4) 25.3 (23.1-27.5) 26.5 (24.2-28.8)"
 [5] "Men 24.9 (23.5-26.2) 24.8 (23.4-26.2) 25.0 (23.6-26.4) 26.8 (25.7-27.8) 26.7 (25.7-27.8) 26.8 (25.7-27.9)"
 [6] "Men 20.6 (19.2-22.0) 20.3 (18.8-21.7) 21.4 (20.0-22.9) 22.6 (21.9-23.4) 21.9 (20.9-22.8) 23.2 (22.3-24.0)"
 [7] "Men 23.3 (22.5-24.0) 22.2 (21.4-23.0) 23.7 (22.9-24.5) 26.2 (25.7-26.7) 25.0 (24.3-25.6) 26.4 (25.9-26.9)"
 [8] "Men 24.0 (22.5-25.4) 23.7 (22.1-25.1) 24.2 (22.7-25.7) 27.1 (26.4-27.7) 26.6 (26.0-27.3) 27.2 (26.5-27.9)"
 [9] "Men 25.0 (23.8-26.4) 24.9 (23.5-26.4) 25.1 (23.8-26.5) 27.0 (25.8-28.1) 27.0 (25.7-28.3) 27.0 (25.7-28.2)"
[10] "Men 20.2 (18.8-21.5) 20.0 (18.6-21.4) 21.4 (19.9-22.9) 22.2 (21.4-23.0) 21.8 (21.0-22.6) 23.1 (22.1-24.1)"
[11] "Men 20.2 (17.7-22.6) 20.1 (17.6-22.5) 21.2 (18.7-23.6) 22.1 (19.8-24.5) 22.0 (19.6-24.3) 23.2 (20.8-25.6)"

Great!

Now for the Women data object

class(Women)

[1] "glue"      "character"

# the female data is already character type
Women<-stringr::str_replace_all(string = Women, 
                                pattern = "na\\*", 
                                replacement = "na\\* na\\*") 
Women[20:30]

 [1] "Women 20.6 (19.1-22.0) 20.1 (18.5-21.6) 21.8 (20.3-23.3) 24.3 (23.5-25.0) 23.2 (22.5-24.0) 25.5 (24.7-26.3)"
 [2] "Women 25.4 (22.2-28.6) na* na* 25.4 (22.2-28.6) 28.5 (25.3-31.6) na* na* 28.5 (25.3-31.6)"                  
 [3] "Women 20.7 (18.4-22.9) 20.3 (18.0-22.6) 23.0 (20.8-25.3) 24.6 (23.7-25.5) 23.8 (22.7-24.9) 25.9 (24.7-27.0)"
 [4] "Women 23.8 (22.3-25.3) 23.0 (21.4-24.5) 24.6 (23.1-26.1) 27.9 (26.5-29.3) 27.0 (25.6-28.5) 28.3 (26.8-29.7)"
 [5] "Women 25.5 (23.5-27.4) 25.7 (23.7-27.7) 25.1 (23.1-27.0) 25.7 (24.4-27.1) 26.0 (24.6-27.5) 25.3 (23.9-26.8)"
 [6] "Women 24.2 (22.2-26.2) 23.7 (21.5-25.8) 25.6 (23.5-27.7) 26.2 (25.3-27.1) 25.0 (23.8-26.2) 27.0 (26.0-28.2)"
 [7] "Women 24.0 (23.2-25.0) 23.1 (22.1-24.1) 24.4 (23.5-25.4) 26.9 (26.2-27.5) 26.4 (25.7-27.2) 26.9 (26.3-27.6)"
 [8] "Women 23.6 (21.4-25.6) 22.9 (20.7-25.0) 24.0 (21.8-26.1) 27.4 (26.6-28.1) 27.0 (26.1-27.8) 27.5 (26.7-28.2)"
 [9] "Women 25.2 (23.4-27.1) 25.5 (23.4-27.6) 25.0 (23.1-26.9) 25.6 (24.1-27.1) 26.0 (24.2-27.8) 25.5 (23.9-27.1)"
[10] "Women 19.9 (18.7-21.2) 19.6 (18.4-20.9) 22.0 (20.7-23.3) 22.2 (21.4-23.1) 21.1 (20.2-22.0) 24.7 (23.7-25.7)"
[11] "Women 19.6 (17.5-21.7) 19.5 (17.4-21.6) 21.7 (19.5-23.8) 21.1 (20.3-22.0) 20.7 (19.8-21.6) 24.0 (23.0-24.9)"

Great, now we can split our data by spaces.

Splitting the data

Men <-tibble::as_tibble(str_split(Men, " ", simplify = TRUE))
# Note here we need to use the column of Women that has the data
head(Men)

# A tibble: 6 x 13
  V1    V2    V3     V4    V5     V6    V7    V8    V9    V10   V11   V12  
  <chr> <chr> <chr>  <chr> <chr>  <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 Men   20.2  (17.8… 19.7  (17.2… 22.4  (20.… 22.8  (20.… 22.5  (20.… 23.6 
2 Men   25.2  (23.9… 25.0  (23.7… 25.4  (24.… 27.0  (26.… 26.9  (25.… 27.0 
3 Men   22.1  (20.8… 21.8  (20.5… 22.3  (21.… 25.1  (24.… 24.8  (24.… 25.2 
4 Men   33.7  (32.7… 32.6  (31.7… 34.0  (32.… 34.3  (33.… 34.6  (33.… 34.2 
5 Men   25.0  (22.5… 25.3  (22.8… 25.0  (22.… 26.8  (24.… 26.8  (24.… 26.8 
6 Men   20.5  (17.9… 20.2  (17.6… 21.4  (18.… 22.6  (20.… 22.0  (19.… 23.2 
# … with 1 more variable: V13 <chr>

Women<-as_tibble(str_split(Women, " ", simplify = TRUE))
# Note here we need to use the column of Women that has the data
head(Women)

# A tibble: 6 x 13
  V1    V2    V3     V4    V5     V6    V7    V8    V9    V10   V11   V12  
  <chr> <chr> <chr>  <chr> <chr>  <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 Women 20.6  (18.4… 20.1  (17.8… 23.2  (20.… 24.4  (23.… 23.6  (22.… 26.3 
2 Women 26.0  (24.1… 26.1  (24.1… 25.9  (23.… 26.0  (24.… 26.2  (24.… 25.9 
3 Women 24.0  (22.2… 23.3  (21.4… 24.8  (22.… 27.4  (26.… 27.0  (26.… 27.5 
4 Women 34.3  (33.1… 34.1  (33.0… 34.4  (33.… 35.3  (33.… 35.0  (33.… 35.4 
5 Women 25.2  (22.0… 25.4  (22.2… 25.2  (22.… 25.3  (22.… 25.2  (21.… 25.3 
6 Women 21.3  (18.0… 20.9  (17.6… 22.7  (19.… 24.4  (21.… 23.3  (20.… 25.7 
# … with 1 more variable: V13 <chr>

Women[20:30,] 

# A tibble: 11 x 13
   V1    V2    V3     V4    V5    V6    V7    V8    V9    V10   V11   V12  
   <chr> <chr> <chr>  <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
 1 Women 20.6  (19.1… 20.1  (18.… 21.8  (20.… 24.3  (23.… 23.2  (22.… 25.5 
 2 Women 25.4  (22.2… na*   na*   25.4  (22.… 28.5  (25.… na*   na*   28.5 
 3 Women 20.7  (18.4… 20.3  (18.… 23.0  (20.… 24.6  (23.… 23.8  (22.… 25.9 
 4 Women 23.8  (22.3… 23.0  (21.… 24.6  (23.… 27.9  (26.… 27.0  (25.… 28.3 
 5 Women 25.5  (23.5… 25.7  (23.… 25.1  (23.… 25.7  (24.… 26.0  (24.… 25.3 
 6 Women 24.2  (22.2… 23.7  (21.… 25.6  (23.… 26.2  (25.… 25.0  (23.… 27.0 
 7 Women 24.0  (23.2… 23.1  (22.… 24.4  (23.… 26.9  (26.… 26.4  (25.… 26.9 
 8 Women 23.6  (21.4… 22.9  (20.… 24.0  (21.… 27.4  (26.… 27.0  (26.… 27.5 
 9 Women 25.2  (23.4… 25.5  (23.… 25.0  (23.… 25.6  (24.… 26.0  (24.… 25.5 
10 Women 19.9  (18.7… 19.6  (18.… 22.0  (20.… 22.2  (21.… 21.1  (20.… 24.7 
11 Women 19.6  (17.5… 19.5  (17.… 21.7  (19.… 21.1  (20.… 20.7  (19.… 24.0 
# … with 1 more variable: V13 <chr>

# We can see that our na values look correct

Looks good!

Adding new header

We can see from our pdf and our object called header what the header was like. Let’s also make a new single line header, but let’s wait to add Country:

#header
new_header <-c("Sex","National_BMI_1985", "National_BMI_1985_CI", "Rural_BMI_1985", 
               "Rural_BMI_1985_CI", "Urban_BMI_1985", "Urban_BMI_1985_CI", "National_BMI_2017",
               "National_BMI_2017_CI","Rural_BMI_2017", "Rural_BMI_2017_CI", "Urban_BMI_2017", 
               "Urban_BMI_2017_CI")

Let’s change the names of our columns of our tibbles to this new header for our Men and Women data

names(Women)<-new_header
names(Men)<-new_header

Now we will add our country data to both our Men and Women tibbles using the bind_cols() function:

Women <- dplyr::bind_cols(country, Women)
# This will add the country as a new column to the Women data object on the left 
Men <- dplyr::bind_cols(country, Men) 
# This will add the country as a new column to the Men data object on the left 
head(Women)

# A tibble: 6 x 14
  country_split Sex   National_BMI_19… National_BMI_19… Rural_BMI_1985
  <chr>         <chr> <chr>            <chr>            <chr>         
1 Afghanistan   Women 20.6             (18.4-22.8)      20.1          
2 Albania       Women 26.0             (24.1-27.9)      26.1          
3 Algeria       Women 24.0             (22.2-25.7)      23.3          
4 American Sam… Women 34.3             (33.1-35.6)      34.1          
5 Andorra       Women 25.2             (22.0-28.4)      25.4          
6 Angola        Women 21.3             (18.0-24.6)      20.9          
# … with 9 more variables: Rural_BMI_1985_CI <chr>, Urban_BMI_1985 <chr>,
#   Urban_BMI_1985_CI <chr>, National_BMI_2017 <chr>,
#   National_BMI_2017_CI <chr>, Rural_BMI_2017 <chr>,
#   Rural_BMI_2017_CI <chr>, Urban_BMI_2017 <chr>, Urban_BMI_2017_CI <chr>

Changing variable names

Here we will change the variable name of the country data to country, currently it is called country_split. Here we will also introduce the concept of the assignment pipe. In this case our pipe operator looks like this %<>% . Using this fancier pipe requires another package, called magrittr. The other simpler pipe options from this package are loaded with tidyverse (if you used library(tidyverse) which loads most tidyverse packages), but not this fancier version.

The > portion of the pipe still behaves like a normal pipe, while the < portion of the pipe makes an assignment to whatever the <is pointing to, just like when we use the typical assignment operator <-.

# library(magrittr) 
# We can't use the `%<>%` unless we load the magrittr package
# We have already done this but we include this for illustrative purposes.
# Here we will just use the traditional assignment strategy

Women<-dplyr::rename(Women, Country = country_split) 
# We have renamed the country_split variable to be called Country

# Here we reassign Men using the pipe strategy
Men %<>%rename(Country = country_split) 
# We have renamed the country_split variable to to be called country
# We have also reassigned Men to the data object 
# which has the country variable renamed

Joining the data

Now we can combine our Men and Women data using the full_join() function of dplyr. This will put all the male data first (x) and all the female data second (y).

BMI<-dplyr::full_join(x = Men, y = Women)
# Let's check the size of our BMI data... it should have 400 rows (obs).
str(BMI)

Classes 'tbl_df', 'tbl' and 'data.frame':   400 obs. of  14 variables:
 $ Country             : chr  "Afghanistan" "Albania" "Algeria" "American Samoa" ...
 $ Sex                 : chr  "Men" "Men" "Men" "Men" ...
 $ National_BMI_1985   : chr  "20.2" "25.2" "22.1" "33.7" ...
 $ National_BMI_1985_CI: chr  "(17.8-22.7)" "(23.9-26.5)" "(20.8-23.3)" "(32.7-34.7)" ...
 $ Rural_BMI_1985      : chr  "19.7" "25.0" "21.8" "32.6" ...
 $ Rural_BMI_1985_CI   : chr  "(17.2-22.2)" "(23.7-26.4)" "(20.5-23.1)" "(31.7-33.5)" ...
 $ Urban_BMI_1985      : chr  "22.4" "25.4" "22.3" "34.0" ...
 $ Urban_BMI_1985_CI   : chr  "(20.0-25.0)" "(24.0-26.7)" "(21.0-23.6)" "(32.9-35.1)" ...
 $ National_BMI_2017   : chr  "22.8" "27.0" "25.1" "34.3" ...
 $ National_BMI_2017_CI: chr  "(20.3-25.3)" "(26.0-27.9)" "(24.5-25.7)" "(33.0-35.6)" ...
 $ Rural_BMI_2017      : chr  "22.5" "26.9" "24.8" "34.6" ...
 $ Rural_BMI_2017_CI   : chr  "(20.0-25.0)" "(25.9-27.9)" "(24.1-25.4)" "(33.1-35.9)" ...
 $ Urban_BMI_2017      : chr  "23.6" "27.0" "25.2" "34.2" ...
 $ Urban_BMI_2017_CI   : chr  "(21.0-26.1)" "(26.0-28.0)" "(24.6-25.9)" "(32.9-35.6)" ...

Sorting the data

Now Lets sort the data by country:

BMI <-dplyr::arrange(BMI, Country)
head(BMI)

# A tibble: 6 x 14
  Country Sex   National_BMI_19… National_BMI_19… Rural_BMI_1985
  <chr>   <chr> <chr>            <chr>            <chr>         
1 Afghan… Men   20.2             (17.8-22.7)      19.7          
2 Afghan… Women 20.6             (18.4-22.8)      20.1          
3 Albania Men   25.2             (23.9-26.5)      25.0          
4 Albania Women 26.0             (24.1-27.9)      26.1          
5 Algeria Men   22.1             (20.8-23.3)      21.8          
6 Algeria Women 24.0             (22.2-25.7)      23.3          
# … with 9 more variables: Rural_BMI_1985_CI <chr>, Urban_BMI_1985 <chr>,
#   Urban_BMI_1985_CI <chr>, National_BMI_2017 <chr>,
#   National_BMI_2017_CI <chr>, Rural_BMI_2017 <chr>,
#   Rural_BMI_2017_CI <chr>, Urban_BMI_2017 <chr>, Urban_BMI_2017_CI <chr>

Our data is looking great! Now we might want to make sure that our observations for each variable look the way we want. In other words, if we want to make plots about National BMI in 1985 then we would need our values to be numeric. Looking at our BMI data using str(), we can see the type of data for each of our variables listed just after variable name and the “:” colon. Looks like none of our BMI data is actually numeric, but of the class character. Let’s change that now.

str(BMI)

Classes 'tbl_df', 'tbl' and 'data.frame':   400 obs. of  14 variables:
 $ Country             : chr  "Afghanistan" "Afghanistan" "Albania" "Albania" ...
 $ Sex                 : chr  "Men" "Women" "Men" "Women" ...
 $ National_BMI_1985   : chr  "20.2" "20.6" "25.2" "26.0" ...
 $ National_BMI_1985_CI: chr  "(17.8-22.7)" "(18.4-22.8)" "(23.9-26.5)" "(24.1-27.9)" ...
 $ Rural_BMI_1985      : chr  "19.7" "20.1" "25.0" "26.1" ...
 $ Rural_BMI_1985_CI   : chr  "(17.2-22.2)" "(17.8-22.4)" "(23.7-26.4)" "(24.1-28.1)" ...
 $ Urban_BMI_1985      : chr  "22.4" "23.2" "25.4" "25.9" ...
 $ Urban_BMI_1985_CI   : chr  "(20.0-25.0)" "(20.9-25.4)" "(24.0-26.7)" "(23.9-27.8)" ...
 $ National_BMI_2017   : chr  "22.8" "24.4" "27.0" "26.0" ...
 $ National_BMI_2017_CI: chr  "(20.3-25.3)" "(23.3-25.4)" "(26.0-27.9)" "(24.8-27.2)" ...
 $ Rural_BMI_2017      : chr  "22.5" "23.6" "26.9" "26.2" ...
 $ Rural_BMI_2017_CI   : chr  "(20.0-25.0)" "(22.5-24.8)" "(25.9-27.9)" "(24.8-27.5)" ...
 $ Urban_BMI_2017      : chr  "23.6" "26.3" "27.0" "25.9" ...
 $ Urban_BMI_2017_CI   : chr  "(21.0-26.1)" "(25.1-27.4)" "(26.0-28.0)" "(24.6-27.2)" ...

We could change these values to be numeric with 6 lines of code like this:

BMI$National_BMI_1985%<>% as.numeric
BMI$Rural_BMI_1985%<>% as.numeric
BMI$Urban_BMI_1985%<>% as.numeric

BMI$National_BMI_2017%<>% as.numeric
BMI$Rural_BMI_2017%<>% as.numeric
BMI$Urban_BMI_2017%<>% as.numeric

#  And if we did this we would see these variables 
#  show as "num" which stands for numeric
str(BMI)

Classes 'tbl_df', 'tbl' and 'data.frame':   400 obs. of  14 variables:
 $ Country             : chr  "Afghanistan" "Afghanistan" "Albania" "Albania" ...
 $ Sex                 : chr  "Men" "Women" "Men" "Women" ...
 $ National_BMI_1985   : num  20.2 20.6 25.2 26 22.1 24 33.7 34.3 25 25.2 ...
 $ National_BMI_1985_CI: chr  "(17.8-22.7)" "(18.4-22.8)" "(23.9-26.5)" "(24.1-27.9)" ...
 $ Rural_BMI_1985      : num  19.7 20.1 25 26.1 21.8 23.3 32.6 34.1 25.3 25.4 ...
 $ Rural_BMI_1985_CI   : chr  "(17.2-22.2)" "(17.8-22.4)" "(23.7-26.4)" "(24.1-28.1)" ...
 $ Urban_BMI_1985      : num  22.4 23.2 25.4 25.9 22.3 24.8 34 34.4 25 25.2 ...
 $ Urban_BMI_1985_CI   : chr  "(20.0-25.0)" "(20.9-25.4)" "(24.0-26.7)" "(23.9-27.8)" ...
 $ National_BMI_2017   : num  22.8 24.4 27 26 25.1 27.4 34.3 35.3 26.8 25.3 ...
 $ National_BMI_2017_CI: chr  "(20.3-25.3)" "(23.3-25.4)" "(26.0-27.9)" "(24.8-27.2)" ...
 $ Rural_BMI_2017      : num  22.5 23.6 26.9 26.2 24.8 27 34.6 35 26.8 25.2 ...
 $ Rural_BMI_2017_CI   : chr  "(20.0-25.0)" "(22.5-24.8)" "(25.9-27.9)" "(24.8-27.5)" ...
 $ Urban_BMI_2017      : num  23.6 26.3 27 25.9 25.2 27.5 34.2 35.4 26.8 25.3 ...
 $ Urban_BMI_2017_CI   : chr  "(21.0-26.1)" "(25.1-27.4)" "(26.0-28.0)" "(24.6-27.2)" ...

Or we can use a more automated way with the map() function of the purrr package:

BMI_numeric <- as_tibble(purrr::map(
  (BMI %>% select(-matches("CI|Sex|Country"))), as.numeric)) 
# The map function of the purrr package allows us to apply the 
# as.numeric() function to all the selected columns of BMI 
# We have selected those that dont contain CI, Sex, or Country in the column name 
# If you are familiar with apply(), this function is very similar
BMI_numeric %<>% mutate(Country =BMI$Country, Sex =BMI$Sex)
head(BMI_numeric)

# A tibble: 6 x 8
  National_BMI_19… Rural_BMI_1985 Urban_BMI_1985 National_BMI_20…
             <dbl>          <dbl>          <dbl>            <dbl>
1             20.2           19.7           22.4             22.8
2             20.6           20.1           23.2             24.4
3             25.2           25             25.4             27  
4             26             26.1           25.9             26  
5             22.1           21.8           22.3             25.1
6             24             23.3           24.8             27.4
# … with 4 more variables: Rural_BMI_2017 <dbl>, Urban_BMI_2017 <dbl>,
#   Country <chr>, Sex <chr>

# In this case our numeric data is class dbl
# This means double-precision floating point numbers (with decimals)
# The alternative numeric option is the integar class

It is generally useful to get the data in what is called long format for other analyses, and particularly for plotting.

For a more detailed description about this, please see this case study

Basically what the long format will have more rows and fewer columns. Thus we can combine or “gather” some columns together. In our case we can put all the different BMI data together in one long column. We will create a new column that tells us what each row of the new BMI column represents. In other words, it will tell us what the original column was.

To do this we will use the gather() function of the tidyr package:

BMI_long <-tidyr::gather(data =BMI_numeric, key = class_BMI, value = BMI, National_BMI_1985:Urban_BMI_2017, factor_key = FALSE) 
# The data value indicates what data you start with
# The key indicates the new name of the column that will now contain 
# information about the values of the columns we will put together
# The value is the name of the column for the values for 
# all the columns that are being put together
# The ":" can be used to indentify the start and end of the column names
# that you are putting together
# The factor_key paramater determines if the newly created key column 
# should be evaluated as a factor or not
head(BMI_long)

# A tibble: 6 x 4
  Country     Sex   class_BMI           BMI
  <chr>       <chr> <chr>             <dbl>
1 Afghanistan Men   National_BMI_1985  20.2
2 Afghanistan Women National_BMI_1985  20.6
3 Albania     Men   National_BMI_1985  25.2
4 Albania     Women National_BMI_1985  26  
5 Algeria     Men   National_BMI_1985  22.1
6 Algeria     Women National_BMI_1985  24  

It would be useful to parse the 1985 and the 2017 data. We can do so by separating the parts of the class_BMI column using the separate() function of the tidyr package. In our case we will replace the class_BMI column with two new columns.

BMI_long%<>% tidyr::separate(class_BMI, c("Region", NA, "Year")) 
# We will separate the class_BMI column of the BMI_long tibble
# This is based on the underscores and create two new columns
# The first column will be called Region
# It will contain the 1st part of the class_BMI data before the 1st underscore
# The 2nd column will be called Year
# It will contain the 3rd part of the data based on underscores
# The middle part of the column will not be used to create any new columns
head(BMI_long)

# A tibble: 6 x 5
  Country     Sex   Region   Year    BMI
  <chr>       <chr> <chr>    <chr> <dbl>
1 Afghanistan Men   National 1985   20.2
2 Afghanistan Women National 1985   20.6
3 Albania     Men   National 1985   25.2
4 Albania     Women National 1985   26  
5 Algeria     Men   National 1985   22.1
6 Algeria     Women National 1985   24  

# Let's see how the dimensions of the BMI data have changed in long format:

dim(BMI_numeric)

[1] 400   8

dim(BMI_long)

[1] 2400    5

Question opportunity: Why exactly are there 2,400 rows now?

Great! our data is very usable now in this format!

Data Exploration

Now that our data is clean and in a format that we can work with, we can start to take a look at the data and how different groups might compare.

Statistically speaking there are tests that can allow us to evaluate if the means of the groups are different. First let’s see how our data looks in general.

General summary

Start here if not working through data import and wrangling

summary(BMI_numeric)

 National_BMI_1985 Rural_BMI_1985 Urban_BMI_1985  National_BMI_2017
 Min.   :18.20     Min.   :17.7   Min.   :19.30   Min.   :20.20    
 1st Qu.:21.50     1st Qu.:20.9   1st Qu.:22.43   1st Qu.:24.07    
 Median :24.00     Median :23.4   Median :24.30   Median :26.15    
 Mean   :23.68     Mean   :23.3   Mean   :24.24   Mean   :26.00    
 3rd Qu.:25.23     3rd Qu.:25.1   3rd Qu.:25.40   3rd Qu.:27.43    
 Max.   :34.30     Max.   :34.1   Max.   :34.40   Max.   :35.30    
                   NA's   :6      NA's   :2                        
 Rural_BMI_2017  Urban_BMI_2017    Country              Sex           
 Min.   :19.80   Min.   :21.50   Length:400         Length:400        
 1st Qu.:23.30   1st Qu.:24.80   Class :character   Class :character  
 Median :26.00   Median :26.30   Mode  :character   Mode  :character  
 Mean   :25.61   Mean   :26.37                                        
 3rd Qu.:27.30   3rd Qu.:27.60                                        
 Max.   :35.00   Max.   :35.40                                        
 NA's   :8       NA's   :2                                            

# We can see that the mean  BMI estimates are larger in 2017
# Is this a meaningful increase?

OK, so our means are fairly similar, but we can see some differences. The question is if those differences are meaningful. This is difficult to determine without a statistical test.

Lets look at the data separately by gender:

#Women:
BMI_numeric %>% 
  subset(Sex == "Women") %>% 
  summary()

 National_BMI_1985 Rural_BMI_1985  Urban_BMI_1985  National_BMI_2017
 Min.   :18.20     Min.   :17.70   Min.   :19.70   Min.   :21.00    
 1st Qu.:22.00     1st Qu.:21.20   1st Qu.:22.85   1st Qu.:24.40    
 Median :24.25     Median :23.90   Median :24.60   Median :26.10    
 Mean   :24.00     Mean   :23.59   Mean   :24.63   Mean   :26.38    
 3rd Qu.:25.50     3rd Qu.:25.40   3rd Qu.:25.60   3rd Qu.:28.00    
 Max.   :34.30     Max.   :34.10   Max.   :34.40   Max.   :35.30    
                   NA's   :3       NA's   :1                        
 Rural_BMI_2017  Urban_BMI_2017    Country              Sex           
 Min.   :20.40   Min.   :21.60   Length:200         Length:200        
 1st Qu.:23.70   1st Qu.:25.30   Class :character   Class :character  
 Median :26.20   Median :26.40   Mode  :character   Mode  :character  
 Mean   :25.98   Mean   :26.84                                        
 3rd Qu.:27.60   3rd Qu.:28.25                                        
 Max.   :35.00   Max.   :35.40                                        
 NA's   :4       NA's   :1                                            

#Men
BMI_numeric %>% 
  subset(Sex == "Men") %>% 
  summary()

 National_BMI_1985 Rural_BMI_1985  Urban_BMI_1985  National_BMI_2017
 Min.   :18.50     Min.   :18.40   Min.   :19.30   Min.   :20.20    
 1st Qu.:21.10     1st Qu.:20.80   1st Qu.:22.10   1st Qu.:23.30    
 Median :23.60     Median :23.20   Median :24.10   Median :26.20    
 Mean   :23.36     Mean   :23.01   Mean   :23.84   Mean   :25.61    
 3rd Qu.:25.00     3rd Qu.:24.90   3rd Qu.:25.10   3rd Qu.:27.10    
 Max.   :33.70     Max.   :32.60   Max.   :34.00   Max.   :34.30    
                   NA's   :3       NA's   :1                        
 Rural_BMI_2017  Urban_BMI_2017   Country              Sex           
 Min.   :19.80   Min.   :21.5   Length:200         Length:200        
 1st Qu.:22.60   1st Qu.:23.8   Class :character   Class :character  
 Median :25.75   Median :26.3   Mode  :character   Mode  :character  
 Mean   :25.24   Mean   :25.9                                        
 3rd Qu.:27.00   3rd Qu.:27.2                                        
 Max.   :34.60   Max.   :34.2                                        
 NA's   :4       NA's   :1                                           

It looks like mean BMIs have increased in all regions for both men and women. It is unclear though if this change is statistically significant.

Distributions

In order to apply a statistical test to compare the means, one of the first things to do is to plot the frequency of the different possible values. This is called a distribution. See here for more informations about distributions. To plot the distribution of a dataset we can use the hist() function to create a histogram.

If the data were what we call normally distributed or followed a Gaussian distribution then the distribution should be equally centered around the mean and would look something like this:

norm_BMI_ex_data <- tibble(norm_data =rnorm(n = 200, mean = 24))
# The rnorm() function allows us to make a vector of normally distributed data
# centered around the mean of 24
head(norm_BMI_ex_data)

# A tibble: 6 x 1
  norm_data
      <dbl>
1      24.0
2      24.2
3      24.6
4      24.1
5      24.4
6      24.9

hist(norm_BMI_ex_data$norm_data)

Alternatively we can plot this using the geom_hist() function of the ggplot2 package. The ggplot2 package creates plots by using layers. Notice in the following code how there is a plus sign between the ggplot() function and the geom_histogram() function. With ggplot2 we select what data we would like to plot using the first function (ggplot()) and then we add on additional layers of complexity (these layers can even involve different data). In this case the geom_histrogram() function initiates the plotting of a histogram with the data that was identified in the ggplot() function. We will see later how we can add many layers to plots with ggplot2. For additional information on using ggplot2, see this case study.

norm_BMI_ex_data %>% 
  ggplot2::ggplot(aes(x = norm_data)) + 
  ggplot2::geom_histogram() 

# We are using the column called norm_data within norm_BMI_ex_data
# With ggplot2 we must first specify the data with ggplot() 
# The geom_histogram() function specifies what type of plot to create

Let’s see how our data looks:

BMI_long%>% 
  ggplot(aes(x = BMI)) + 
  geom_histogram() 

OK, so the data looks like it is what we call right-skewed because the tail of the distribution is longer on the right side (in other words, it is a bit wider on that side), but it looks fairly normally distributed. See here for more information about skewed distributions.

If we parse our data, what does it look like? The easiest way to do this is to use some other functions of the ggplot2 package, particularly the facet_wrap() function, which will allow us to look at differences in the distribution of the BMI data by year, gender, and region. We can sequentially divide our plots by deeper levels using multiple variables and the + plus sign within facet_wrap().

BMI_long %>%
  ggplot(aes(x=BMI)) +
  geom_histogram() +
  facet_wrap(~ Sex)

BMI_long %>%
  ggplot(aes(x=BMI)) +
  geom_histogram() +
  facet_wrap(~ Sex + Year) 

BMI_long %>%
  ggplot(aes(x=BMI)) +
  geom_histogram() +
  facet_wrap(~ Sex + Year + Region) 

OK, so some of these plots look pretty similar to our normal distribution plot, like the Women 2017 Urban plot (although it is right-skewed). However some of the other plots, like the Men 2017 Urban plot, look very different. These plots show what is called a bimodal distribution, meaning that there appear to be two peaks, or in other words two different values (and values near each of these) that are the most frequent.

Quantile-Quantile Plots

Let’s use a method called the quantile-quantile plot to determine if the data is indeed normally distributed. These plots are called Q-Q plots. This method allows us to test the fit of known theoretical distributions (like the normal distribution) with our observed distribution. To do this we will plot the quantiles of our data on the y-axis and the quantiles of the theoretical normal distribution on the x-axis. If the quantiles line up then we can say that our data is fairly normal. What exactly is a quantile? This is a division of the data distribution into roughly equal portions. See here for more information about Q-Q plots.

Here is an example of a Q-Q plot for the normally distributed data that we just created using the stat_qq() function of the ggplot2 package:

norm_BMI_ex_data %>%
  ggplot(aes(sample =norm_data)) +
  stat_qq() + # to add the line we need to also include stat_qq_line()
  stat_qq_line()

Here we can see that the quantiles are fairly similar between the observed and theoretical data. We see that the points mostly fall on the line, however there are some points that are a bit further from the line as we get to the extreme quantiles. Notice that the sample quantiles (which will be fairly similar to our real BMI data quantiles) on the y-axis has the same range as the values that we created. So values that are bellow 22 for example are represented as the points bellow 22 on the y-axis. As expected we see that about half the points are bellow our mean of 24.

If we were to use different data that had a range of different values our y-axis would shift according to the range of values. For example the Orange data within the installation of R includes data about the circumference of orange trees in millimeters. Here we can see that the quantiles are quite different but reflect the range of orange tree circumferences.

range(Orange$circumference)

[1]  30 214

Orange%>%
  ggplot(aes(sample =circumference)) +
  stat_qq() + # to add the line we need to also include stat_qq_line()
  stat_qq_line()

Let’s take a look at our BMI data:

BMI_long %>%
  ggplot(aes(sample =BMI)) +
  stat_qq() + # to add the line we need to also include stat_qq_line()
  stat_qq_line() +
  facet_wrap(~ Sex + Year + Region)

This is a little hard to see, some plots look pretty good, but while others seem to be skewed, lets take a closer look.

For the sake of simplicity, we are going to focus on the data from the women. This is because women were identified in the paper to have larger increases in BMI. If we want to perform tests on these groups (women in each region and year) specifically, then we need to look at the distribution of each of these subsets.

BMI_long %>%
  filter(Sex == "Women") %>%
  ggplot(aes(sample =BMI)) +
  stat_qq() + # to add the line we need to also include stat_qq_line()
  stat_qq_line() +
  facet_wrap(~ Year + Region)

We can see that at the extremes of our quantiles, for most of our data, our tails are not very similar to the theoretical distribution. We can can also determine that our data is right- or positive-skewed because the extreme values are above the line that indicates a perfect match between the theoretical and observed data.

The rural data looks more normal than the urban data, but if we were to use statistical tests that rely on normality, this would require that both groups are normally distributed.

Finally, some statisticians also use the Shapiro-Wilk test for normality when the Q-Q plot is a bit unclear. Many statisticians however would conclude from the Q-Q plots that normality appears to be violated in our data.

Shapiro-Wilk test

For illustrative purposes we will show how we can use the dplyr summarize() and group_by() functions to perform the Shapiro-Wilk test for all the subsets we are interested in.

norm_BMI_ex_data$norm_data %>% 
  shapiro.test() 

    Shapiro-Wilk normality test

data:  .
W = 0.99174, p-value = 0.3151

# example of the test output for normal data - the p value is not < 0.05

BMI_long %>% 
  filter(Sex == "Women") %>% 
  group_by(Year, Region) %>% 
  summarize(shapiro_test = shapiro.test(BMI)$p.value)

# A tibble: 6 x 3
# Groups:   Year [2]
  Year  Region   shapiro_test
  <chr> <chr>           <dbl>
1 1985  National  0.0000478  
2 1985  Rural     0.000679   
3 1985  Urban     0.000000332
4 2017  National  0.00363    
5 2017  Rural     0.0108     
6 2017  Urban     0.0000130  

We see that all the data does not appear to be normally distributed.

Data Analysis

We can use statistical tests like the student t-test to determine if two means are significantly different.

With the t-test, we perform what is called two means hypothesis tests . We define what is called the null hypothesis , which is what we assume by default or the baseline, that there is no difference in the two means:

                                    Ho: μ1 = μ2
                                    
                           μ1 is the true mean of one group
                           μ2 is the true mean of the other group
                           

Thus like in law, we accept this null hypothesis, unless we have enough evidence to suggest that we should reject it - similarly to the idea that we by default assume that indivudals are not guilty until proven otherwise. If we reject the null hypothesis, we then accept an alternative hypothesis , this can vary for different statistical tests, but in the case of the t-test we evaluate if the two means are not equal:

                                    Ha: μ1 ≠ μ2
                                    

Remeber that we can’t know the true mean for our populations of interest, instead we approximate or estimate our means based on the sample that we have in our data.

We are interested in comparing the means of female rural and urban BMI measurements for both years.

There are two possible classes of statistical tests that we could run to compare the means of these two groups:

1. Parametric
2. Nonparametric

Parametric tests are based on assumptions about the distribution of the data, while Nonparametric tests do not rely on this assumption. It is called “parametric” because aspects about the distribution of the data like the mean aer called parameters when we describe a population. In parametric tests we estimate the parameters of the true population of interest using a sample of that population. These estiamtes are called statistics. See here for more information about the difference between these two classes of tests.

Parametric two sample mean tests

Often when comparing two groups we might perform a two sample t-test to determine if the means of each group is different. The two sample t-test however, relies on several assumptions:

1. The data for both groups is normally distributed
2. The variance of both groups is similar
3. The number of observations is similar for both groups - thus they are balanced
4. The observations are independent (meaning that observations do not influence eachother)

If these assumptions are violated, this doesn’t necessarily mean we can’t perform a t-test. It just means we may need to consider the following options:

1. Transformation of the data to make it more normally distributed
2. Welch’s t-test also call the unequal variance t-test we may need to modify the way we perform the t-test to account for the difference in the variance in the two groups
3. Permutation/resampling methods to deal with violations of normality or imbalance.

Alternatively, we can use a nonparametric test like the Wilcoxon–Mann–Whitney (WMW) test. These tests are often a good option when multiple assumptions are violated are when sample sizes are small. We will explore these options.

Our data has a balance of observations for both groups - in fact they are equal, thus that assumption is not violated. If it were violated, we would want to consider using permutation methods which are also a good option for violations of normality. To learn more about these methods see here.

If we needed to check if our samples were imbalanced, we could use the table() function:

table(BMI_long$Sex, BMI_long$Region, BMI_long$Year)

, ,  = 1985

       
        National Rural Urban
  Men        200   200   200
  Women      200   200   200

, ,  = 2017

       
        National Rural Urban
  Men        200   200   200
  Women      200   200   200

We can see that the number of observations for each possible group of interest is the same.

The t-test is also fairly robust to non-normality if the data is relatively large, due to what is called the central limit theorem.

We have an n of 200, which should be sufficient but let’s investigate the nonparametric tests further.

Often we would check if the variance of the rural and urban data is equal using the var.test() function. However this is an F test and assumes that the data is normally distributed. Instead we will use the mood.test() function which performs the Mood’s two-sample test for a difference in scale parameters and does not assume that the data is normally distributed. We will also introduce the pull() function of the dplyr package.

#pull() will allow us to take the BMI column after we filter BMI_long

mood.test(dplyr::pull(filter(BMI_long, 
                             Sex == "Women", 
                             Year == "2017", 
                             Region == "Rural"), BMI), 
          dplyr::pull(filter(BMI_long,
                             Sex == "Women", 
                             Year == "2017", 
                             Region == "Urban"), BMI))

    Mood two-sample test of scale

data:  dplyr::pull(filter(BMI_long, Sex == "Women", Year == "2017",  and dplyr::pull(filter(BMI_long, Sex == "Women", Year == "2017",     Region == "Rural"), BMI) and     Region == "Urban"), BMI)
Z = 2.9189, p-value = 0.003513
alternative hypothesis: two.sided

# p value <.05, conclude that variance is not equal
# reject the null: no difference in the variance of the distributions

mood.test(pull(filter(BMI_long, Sex == "Women", 
                      Year == "1985", 
                      Region == "Rural"), BMI), 
          pull(filter(BMI_long, 
                      Sex == "Women", 
                      Year == "1985", 
                      Region == "Urban"), BMI))

    Mood two-sample test of scale

data:  pull(filter(BMI_long, Sex == "Women", Year == "1985", Region ==  and pull(filter(BMI_long, Sex == "Women", Year == "1985", Region ==     "Rural"), BMI) and     "Urban"), BMI)
Z = 3.1305, p-value = 0.001745
alternative hypothesis: two.sided

# p value <.05, conclude that variance is not equal
# reject the null: no difference in the variance of the distributions

mood.test(pull(filter(BMI_long, 
                      Sex == "Women", 
                      Year == "1985", 
                      Region == "Rural"), BMI), 
          pull(filter(BMI_long, Sex == "Women", 
                      Year == "2017", 
                      Region == "Rural"), BMI))

    Mood two-sample test of scale

data:  pull(filter(BMI_long, Sex == "Women", Year == "1985", Region ==  and pull(filter(BMI_long, Sex == "Women", Year == "2017", Region ==     "Rural"), BMI) and     "Rural"), BMI)
Z = -0.24228, p-value = 0.8086
alternative hypothesis: two.sided

# p value >.05, conclude that variance is equal
# accept the null: no difference in the variance of the distributions

mood.test(pull(filter(BMI_long, 
                      Sex == "Women", 
                      Year == "1985", 
                      Region == "Urban"), BMI), 
          pull(filter(BMI_long, Sex == "Women", 
                      Year == "2017", 
                      Region == "Urban"), BMI))

    Mood two-sample test of scale

data:  pull(filter(BMI_long, Sex == "Women", Year == "1985", Region ==  and pull(filter(BMI_long, Sex == "Women", Year == "2017", Region ==     "Urban"), BMI) and     "Urban"), BMI)
Z = 1.5317, p-value = 0.1256
alternative hypothesis: two.sided

# p value >.05, conclude that variance is equal
# accept the null: no difference in the variance of the distributions

Our p value is less than .05 for both tests, thus we reject our null hypothesis that there is no difference in the variance. Therefore, we conclude that the variance is not equal and that our data also violates this assumption.

We will perform a special t.test where we account for the fact that our variance is not equal.

Another important consideration is that the data is what we call paired, meaning that the measurements from the rural and urban areas are not independent. That is because we have a rural and urban measurement mean for nearly every country. Thus these values may be more similar to one another if they come from the same country. This is also true for the male and female measurements from the same country or the values in the same countries from 1985 and later in 2017. However, we are assuming that measurements between different countries are independent, thus this assumption is not violated, making it reasonable to peform the welch’s or paired t-test.

When we perform a paired t-test our hypothesis is slightly different from the typical student’s t-test. In this case we are testing the differences among the pairs of observations and how close these differences are to zero.In this case our null hypothesis is that the mean differences is equal to zero:

                                    Ho: μd = 0
                                    
                           μd is the true mean differences 
                           between paired observations of the two groups

In this case the alternative hypthesis is that the mean of the differences is not equal to zero

                                    Ha: μd ≠ 0
                                    
                           μd is the true mean differences 
                           between paired observations of the two groups

t.test(pull(filter(BMI_long, Sex == "Women", 
                   Year == "2017", 
                   Region == "Rural"), BMI), 
       pull(filter(BMI_long, Sex == "Women", 
                   Year == "2017", 
                   Region == "Urban"), BMI), 
       var.equal = FALSE, paired = TRUE)

    Paired t-test

data:  pull(filter(BMI_long, Sex == "Women", Year == "2017", Region ==  and pull(filter(BMI_long, Sex == "Women", Year == "2017", Region ==     "Rural"), BMI) and     "Urban"), BMI)
t = -10.356, df = 194, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.0573625 -0.7190478
sample estimates:
mean of the differences 
             -0.8882051 

# means are different - p value <.05 reject the null: no difference in the means

t.test(pull(filter(BMI_long, Sex == "Women", 
                   Year == "1985",
                   Region == "Rural"), BMI), 
       pull(filter(BMI_long, Sex == "Women", 
                   Year == "1985", 
                   Region == "Urban"), BMI),
       var.equal = FALSE, paired = TRUE)

    Paired t-test

data:  pull(filter(BMI_long, Sex == "Women", Year == "1985", Region ==  and pull(filter(BMI_long, Sex == "Women", Year == "1985", Region ==     "Rural"), BMI) and     "Urban"), BMI)
t = -14.095, df = 195, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.1870263 -0.8956268
sample estimates:
mean of the differences 
              -1.041327 

# means are different - p value <.05 reject the null: no difference in the means

t.test(pull(filter(BMI_long, Sex == "Women", 
                   Year == "1985", 
                   Region == "Rural"), BMI), 
       pull(filter(BMI_long, Sex == "Women", 
                   Year == "2017", 
                   Region == "Rural"), BMI),
       var.equal = TRUE, paired = TRUE)

    Paired t-test

data:  pull(filter(BMI_long, Sex == "Women", Year == "1985", Region ==  and pull(filter(BMI_long, Sex == "Women", Year == "2017", Region ==     "Rural"), BMI) and     "Rural"), BMI)
t = -22.119, df = 195, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.591762 -2.167422
sample estimates:
mean of the differences 
              -2.379592 

# means are different - p value <.05 reject the null: no difference in the means

t.test(pull(filter(BMI_long, Sex == "Women", 
                   Year == "1985", 
                   Region == "Urban"), BMI), 
       pull(filter(BMI_long, Sex == "Women", 
                   Year == "2017", 
                   Region == "Urban"), BMI),
       var.equal = TRUE, paired = TRUE)

    Paired t-test

data:  pull(filter(BMI_long, Sex == "Women", Year == "1985", Region ==  and pull(filter(BMI_long, Sex == "Women", Year == "2017", Region ==     "Urban"), BMI) and     "Urban"), BMI)
t = -24.378, df = 198, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.383938 -2.027118
sample estimates:
mean of the differences 
              -2.205528 

# means are different - p value <.05 reject the null: no difference in the means
# the t value shows us that the directionality of the difference
# the t value was negative thus BMI was lower in 1985 (1st group) 
# compared to 2017 (2nd group)

Question opportunity: Looking at the t value, was global BMI lower in Rural or Urban areas in 1985?

Now we will try transform our data to make it more normally distributed. One way to do this is to take the logarithm of the data values. Then we will see how this influences the results. Again we will focus on the data for women.

# create a new column with the log version of the BMI variable using mutate()

BMI_long_log <- mutate(BMI_long, log_BMI=log(pull(BMI_long,BMI)))

BMI_long_log %>%
  filter(Sex == "Women") %>%
  ggplot(aes(sample =log_BMI)) +
  stat_qq() + # to add the line we need to also include stat_qq_line()
  stat_qq_line() +
  facet_wrap(~ Year + Region)

BMI_long_log %>% 
  filter(Sex == "Women") %>% 
  group_by(Year, Region) %>% 
  summarize(shapiro_test = shapiro.test(log_BMI)$p.value)

# A tibble: 6 x 3
# Groups:   Year [2]
  Year  Region   shapiro_test
  <chr> <chr>           <dbl>
1 1985  National     0.00315 
2 1985  Rural        0.0105  
3 1985  Urban        0.000334
4 2017  National     0.293   
5 2017  Rural        0.0784  
6 2017  Urban        0.00416 

# recall what it was before
BMI_long %>% 
  filter(Sex == "Women") %>% 
  group_by(Year, Region) %>% 
  summarize(shapiro_test = shapiro.test(BMI)$p.value)

# A tibble: 6 x 3
# Groups:   Year [2]
  Year  Region   shapiro_test
  <chr> <chr>           <dbl>
1 1985  National  0.0000478  
2 1985  Rural     0.000679   
3 1985  Urban     0.000000332
4 2017  National  0.00363    
5 2017  Rural     0.0108     
6 2017  Urban     0.0000130  

The data appears to be more similar to the normal distribution, although not quite. Again, our sample size of 200 is quite large and the t-test is generally quite robust to violations of normality with large n, thus the modified t-test to account for unequal variance might be a good option using the log normalized data, as it is at least more normally distributed.

Let’s see the results of the t-test with the transformed data:

t.test(pull(filter(BMI_long_log, Sex == "Women", 
                   Year == "2017", 
                   Region == "Rural"), log_BMI), 
       pull(filter(BMI_long_log, Sex == "Women", 
                   Year == "2017", 
                   Region == "Urban"), log_BMI), 
       var.equal = FALSE, paired = TRUE)

    Paired t-test

data:  pull(filter(BMI_long_log, Sex == "Women", Year == "2017", Region ==  and pull(filter(BMI_long_log, Sex == "Women", Year == "2017", Region ==     "Rural"), log_BMI) and     "Urban"), log_BMI)
t = -10.058, df = 194, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.04242774 -0.02851589
sample estimates:
mean of the differences 
            -0.03547182 

# means are different - p value <.05 reject the null: no difference in the means

t.test(pull(filter(BMI_long_log, Sex == "Women", 
                   Year == "1985", 
                   Region == "Rural"), log_BMI), 
       pull(filter(BMI_long_log, Sex == "Women", 
                   Year == "1985", 
                   Region == "Urban"), log_BMI),
       var.equal = FALSE, paired = TRUE)

    Paired t-test

data:  pull(filter(BMI_long_log, Sex == "Women", Year == "1985", Region ==  and pull(filter(BMI_long_log, Sex == "Women", Year == "1985", Region ==     "Rural"), log_BMI) and     "Urban"), log_BMI)
t = -13.962, df = 195, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.05214677 -0.03923811
sample estimates:
mean of the differences 
            -0.04569244 

# means are different - p value <.05 reject the null: no difference in the means

t.test(pull(filter(BMI_long_log, Sex == "Women", 
                   Year == "1985", 
                   Region == "Rural"), log_BMI), 
       pull(filter(BMI_long_log, Sex == "Women", 
                   Year == "2017", 
                   Region == "Rural"), log_BMI),
       var.equal = TRUE, paired = TRUE)

    Paired t-test

data:  pull(filter(BMI_long_log, Sex == "Women", Year == "1985", Region ==  and pull(filter(BMI_long_log, Sex == "Women", Year == "2017", Region ==     "Rural"), log_BMI) and     "Rural"), log_BMI)
t = -22.369, df = 195, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.10617051 -0.08896626
sample estimates:
mean of the differences 
            -0.09756839 

# means are different - p value <.05 reject the null: no difference in the means

t.test(pull(filter(BMI_long_log, Sex == "Women", 
                   Year == "1985", 
                   Region == "Urban"), log_BMI), 
       pull(filter(BMI_long_log, Sex == "Women", 
                   Year == "2017", 
                   Region == "Urban"), log_BMI),
       var.equal = TRUE, paired = TRUE)

    Paired t-test

data:  pull(filter(BMI_long_log, Sex == "Women", Year == "1985", Region ==  and pull(filter(BMI_long_log, Sex == "Women", Year == "2017", Region ==     "Urban"), log_BMI) and     "Urban"), log_BMI)
t = -23.977, df = 198, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.09377834 -0.07952498
sample estimates:
mean of the differences 
            -0.08665166 

# means are different - p value <.05 reject the null: no difference in the means

We can see that our results are quite similar to that of the original data, however the t values are slightly smaller. In other cases we may see a much more dramatic influence of transforming our data.

Now, let’s take a look at nonparametric tests, which are also a great option when the assumptions of the t-test are violated.

Nonparametric two sample tests

There are two nonparametric options to consider when the assumptions of the t-test are violated. The Wilcoxon signed rank test (for paired data - the alternative is Wilcoxon rank sum test (also called the Mann-Whitney U test) for independent samples) and the two-sample Kolmogorov-Smirnov test both do not assume normality (has both paired and unpaired methods). Thus these tests should be considered when the data of either groups does not appear to be normally distributed and particularly when the number of samples is low.

Importantly the two-sample Kolmogorov-Smirnov (KS) test does not assume normality or equal variance, while the Wilcoxon signed rank test does assume equal variance. Here is how you would perform these tests. However in our case, because the variance is not equal between some of our groups of interest, the KS test would be more appropriate. Both the t-test and the KS test evaluate if the distributions of the two groups are identical, however the KS test does not particularly test any aspect of the distribution like the mean, therefore there are no confidence intervals in the output using this test.

ks.test(pull(filter(BMI_long, Sex == "Women", 
                    Year == "2017", 
                    Region == "Rural"), BMI), 
        pull(filter(BMI_long, Sex == "Women", 
                    Year == "2017", 
                    Region == "Urban"), BMI),
        paired = TRUE)

    Two-sample Kolmogorov-Smirnov test

data:  pull(filter(BMI_long, Sex == "Women", Year == "2017", Region ==  and pull(filter(BMI_long, Sex == "Women", Year == "2017", Region ==     "Rural"), BMI) and     "Urban"), BMI)
D = 0.20006, p-value = 0.0007385
alternative hypothesis: two-sided

ks.test(pull(filter(BMI_long, Sex == "Women", 
                    Year == "1985", 
                    Region == "Rural"), BMI), 
        pull(filter(BMI_long, Sex == "Women", 
                    Year == "1985", 
                    Region == "Urban"), BMI),
        paired = TRUE)

    Two-sample Kolmogorov-Smirnov test

data:  pull(filter(BMI_long, Sex == "Women", Year == "1985", Region ==  and pull(filter(BMI_long, Sex == "Women", Year == "1985", Region ==     "Rural"), BMI) and     "Urban"), BMI)
D = 0.19914, p-value = 0.0007779
alternative hypothesis: two-sided

What about the difference in female BMI from 1985 to 2017 for both regions? Recall that the variance was equal for these comparisons.

wilcox.test(pull(filter(BMI_long, Sex == "Women", 
                        Year == "1985", 
                        Region == "Rural"), BMI),
           pull(filter(BMI_long, Sex == "Women", 
                       Year == "2017", 
                       Region == "Rural"), BMI),
           paired = TRUE)

    Wilcoxon signed rank test with continuity correction

data:  pull(filter(BMI_long, Sex == "Women", Year == "1985", Region ==  and pull(filter(BMI_long, Sex == "Women", Year == "2017", Region ==     "Rural"), BMI) and     "Rural"), BMI)
V = 273, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0

wilcox.test(pull(filter(BMI_long, Sex == "Women", 
                        Year == "1985", 
                        Region == "Urban"), BMI),
           pull(filter(BMI_long, Sex == "Women", 
                       Year == "2017", 
                       Region == "Urban"), BMI),
           paired = TRUE)

    Wilcoxon signed rank test with continuity correction

data:  pull(filter(BMI_long, Sex == "Women", Year == "1985", Region ==  and pull(filter(BMI_long, Sex == "Women", Year == "2017", Region ==     "Urban"), BMI) and     "Urban"), BMI)
V = 189, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0

There is a significant difference across time for both regions, as we saw with the t-test. There is also a significant difference by region for each year. However, the p-values are a bit larger for the KS test results than we saw with the t-test.

Multiple Testing Correction

It’s important to note that ultimately we wanted to test 4 different tests:

1. Is there a difference in mean BMI for women between Rural and Urban areas in 1985?
2. Is there a difference in mean BMI for women between Rural and Urban areas in 2017?
3. Is there a difference in mean BMI for women between 1985 and 2017 in Rural areas?
4. Is there a difference in mean BMI for women between 1985 and 2017 in Urban areas?

When we test multiple questions like this, we need to correct for the multiple tests that we are performing.

The more we test data, the more likely we are to see a significant finding just by chance.

One way to do this is using the Bonferroni method.

In this method we would divide our significance threshold (generally 0.05) by the number of tests.

.05/4 # 4 tests

[1] 0.0125

Our new significance threshold is now 0.0125. Thus our p-values should be less than this value for us to reject the null that there is no difference in means. In all cases, our p-values were less than 0.0125. So we see a significant difference in the means of the groups after multiple testing correction for our different tests.

Again, it would be reasonable to use the t-test because it is robust to deviations in normality when samples are relatively large. We can see that we obtained the same results regardless of the test that we used. However, if sample sizes are small (generally speaking n<15 for each group), then these nonparametric options are useful to know.

The D values in the output of our KS tests, show the magnitude of distance in the difference between the distributions of the groups tested. You may notice that the D value is larger for the tests of BMI across time rather then across region. In these tests the p-value was also smaller.

Data Visualization

Again we will utilize ggplot2 to create plots to look at the directional and magnitude of the differences in BMI we are interested in. If you need additional information please see here. The top two lines of the code for the following plots, filter the data to only specific values of interest. Then we layer what is called a jitter on top of a box plot. A jitter is essentially a dot plot but with some variation on the location of the points so that they do not line up vertically which can make the individual points difficult to see.

Let’s look at the national mean BMI estimates for each of the years:

BMI_long %>% 
  filter(Sex == "Women", 
         Year == "1985", 
         Region != "National") %>%
  ggplot(aes(x = Region, y = BMI)) +
  geom_boxplot() +
  geom_jitter(width = .3) 

# the width determines how wide the points are plotted

BMI_long %>% 
  filter(Sex == "Women", 
         Year == "2017", 
         Region != "National") %>%
  ggplot(aes(x = Region, y = BMI)) +
  geom_boxplot() +
  geom_jitter(width = .3)

Now Let’s look at the change in rural and urban mean BMI estimates:

#this time we will add titles for Rural and Urban using ggtitle()
BMI_long %>% 
  filter(Sex == "Women", 
         Year %in% c("1985", "2017"), 
         Region == "Rural") %>%
  # filtering this way allows us to keep data from both years
  ggplot(aes(x = Year, y = BMI)) +
  geom_boxplot() +
  geom_jitter(width = .3) +
  ggtitle("Rural")

BMI_long %>% 
  filter(Sex == "Women", 
         Year %in% c("1985", "2017"), 
         Region == "Urban") %>%
  ggplot(aes(x = Year, y = BMI)) +
  geom_boxplot() +
  geom_jitter(width = .3) +
  ggtitle("Urban")

Let’s put the plots together to see how the change over the years differs between the regions. We will use again use facet_wrap() to do this:

BMI_long %>% 
  filter(Sex == "Women", 
         Year %in% c("1985", "2017"), 
         Region %in% c("Rural", "Urban")) %>%
  # filtering this way allows us to keep data from both years
  ggplot(aes(x = Year, y = BMI)) +
  geom_boxplot() +
  geom_jitter(width = .3) +
  facet_wrap(~ Region)

# we would also get national data without `Region %in% c("Rural", "Urban"))`

Indeed the change in BMI over time looks bigger in the rural areas!

How do the different countries compare? Or in other words what do the individual dots represent in our box plots? We will take a look using geom_label():

BMI_long %>% 
  filter(Sex == "Women", 
         Year %in% c("1985", "2017"), 
         Region == "Rural") %>%
  ggplot(aes(x = Year, y = BMI)) +
  geom_boxplot() +
  geom_jitter(width = .3) +
  geom_label(aes(label = Country))

If we include all country names this is a bit too much… so perhaps we should focus on just the extreme BMI values using filter() function of thedplyr package. We will also use the ggrepel package to have our labels not overlap each other.

BMI_long %>% 
  filter(Sex == "Women", 
         Year %in% c("1985", "2017"), 
         Region == "Rural") %>%
  ggplot(aes(x = Year, y = BMI)) +
  geom_boxplot() +
  geom_jitter(data=BMI_long %>% 
              filter(Sex == "Women", 
                     Year %in% c("1985", "2017"), 
                     Region == "Rural") %>%
              filter(BMI>31 | BMI<19 | Country == "United States of America"), 
              aes(x =Year, y =BMI),
              width = .02) +
  ggrepel::geom_text_repel(data=BMI_long %>% 
              filter(Sex == "Women", 
                     Year %in% c("1985", "2017"), 
                     Region == "Rural") %>%
              filter(BMI>31 | BMI<19 | Country == "United States of America"), 
              aes(x =Year, y =BMI, label = Country))

And let’s fill the box plots with color and outline in black:

BMI_long %>% 
  filter(Sex == "Women", 
         Year %in% c("1985", "2017"), 
         Region == "Rural") %>%
  ggplot(aes(x = Year, y = BMI)) +
  geom_boxplot(outlier.shape = NA, color = "black", aes(fill = Year)) +
  geom_jitter(data=BMI_long %>% 
              filter(Sex == "Women", 
                     Year %in% c("1985", "2017"), 
                     Region == "Rural") %>%
              subset(BMI>31 | BMI<19 | Country == "United States of America"), 
              aes(x =Year, y =BMI),
              width = .02) +
  ggrepel::geom_text_repel(data=BMI_long %>% 
              filter(Sex == "Women", 
                     Year %in% c("1985", "2017"), 
                     Region == "Rural") %>%
              subset(BMI>31 | BMI<19 | Country == "United States of America"), 
              aes(x =Year, y =BMI, label = Country))

Overall differences

Let’s take a look at all the data together:

ggplot(BMI_long, aes(x = Year, y = BMI, col = Region)) +
  geom_boxplot(outlier.shape = NA, color = "black" , aes(fill = Region)) + 
  facet_grid(~ Sex)   

That’s useful, but let’s look at the individual points and include our country labels, to do so lets change our United States of America label to USA:

BMI_long$Country <-BMI_long$Country %>%
  str_replace( pattern = "United States of America", replacement = "USA")


ggplot(BMI_long, aes(x = Year, y = BMI, col = Year)) + 
  geom_boxplot(outlier.shape = NA, color = "black" , aes(fill = Year)) + 
  facet_grid(~ Sex + Region) + 
  geom_hline(yintercept=30, linetype="dashed", color = "red", size =1) + 
  # This will add a horizontal dashed line to indicate the obesity BMI threshold
  geom_jitter(data=subset(BMI_long), 
              aes(x =Year, y =BMI), 
              width = .2, size =2, shape =21, color = "black", fill = "gray") + 
  # This will add the individual country data points
  # The shape 21 allows for a different fill and outline color
  # The width determines how wide the jitter points are plotted
  geom_jitter(data=subset(BMI_long, Country == "USA"), 
              aes(x =Year, y =BMI), 
              width = .02, size =12, shape =21, color = "black", fill = "gray") + 
  # This will add points that are larger for the USA data
  geom_text(data=subset(BMI_long,  Country == "USA"), 
            aes(x =Year, y =BMI,label=Country), 
            color = "black") + 
  # This will add USA labels to the USA points
  theme(legend.position = "none", 
  # This is useful for removing the legend
        axis.text.x = element_text(size = 15,angle = 30), 
        # this changes the size and angle of the x axis point labels 
        axis.text.y = element_text(size = 20), 
        axis.title.y = element_text(size =15), 
        axis.title.x = element_text(size =15), 
        strip.text.x = element_text(size = 15)) +
  # This changes the size of x axis labels for the facet 
        ggtitle( "Differences in BMI Over Time and Across Region Type and Gender") 

  # Add a  plot title

Here we can see that overall BMI appears to be increasing globally over time. Additionally we can see that this is occurring not just in urban areas, but also in rural areas. The US is consistently above the median in all strata of the data. In general, the female data shows higher BMI values than the male data. The rural USA BMI appears to be higher than the urban BMI for both men and women. Many countries have average BMI estimates above the obesity threshold of 30. Thus it appears that education and outreach programs for weight management should focus on both rural and urban areas and both genders. Education and assistance for women may be especially helpful.

Differences in rate of change in BMI

How does the rate of change in BMI differ between groups? Which group might especially need attention?

First let’s calculate the differences in BMI from 2017 and 1985 and add this to our BMI_long data object:

BMI_numeric$Rural_difference <- 
     BMI_numeric$Rural_BMI_2017 - BMI_numeric$Rural_BMI_1985

BMI_numeric$Urban_difference <- 
     BMI_numeric$Urban_BMI_2017 - BMI_numeric$Urban_BMI_1985

BMI_numeric$National_difference <- 
     BMI_numeric$National_BMI_2017 - BMI_numeric$National_BMI_1985

BMI_diff_long <- BMI_numeric %>% 
  select(Country: National_difference) %>% 
gather(key = Type, value = Difference , Rural_difference:National_difference)

head(BMI_diff_long)

# A tibble: 6 x 4
  Country     Sex   Type             Difference
  <chr>       <chr> <chr>                 <dbl>
1 Afghanistan Men   Rural_difference     2.8   
2 Afghanistan Women Rural_difference     3.5   
3 Albania     Men   Rural_difference     1.90  
4 Albania     Women Rural_difference     0.1000
5 Algeria     Men   Rural_difference     3     
6 Algeria     Women Rural_difference     3.70  

Let’s replace “United states of America” with “USA” and make a plot with this data to compare the change in BMI:

BMI_diff_long$Country <-BMI_diff_long$Country %>%
  str_replace( pattern = "United States of America", replacement = "USA")


ggplot(BMI_diff_long, aes(x = Type, y = Difference, col = Type)) + 
  geom_boxplot(outlier.shape = NA, color = "black" , aes(fill = Type)) + 
  facet_grid(~ Sex) +
  geom_jitter(data=subset(BMI_diff_long), 
              aes(x =Type, y =Difference), 
              width = .2, size =2, shape =21, color = "black", fill = "gray") + 
  # This will add the individual country data points
  # The shape 21 allows for a different fill and outline color
  # The width determines how wide the jitter points are plotted
  geom_jitter(data=subset(BMI_diff_long, Country == "USA"), 
              aes(x =Type, y =Difference), 
              width = .02, size =12, shape =21, color = "black", fill = "gray") + 
  # This will add points that are larger for the USA data
  geom_text(data=subset(BMI_diff_long,  Country == "USA"), 
            aes(x =Type, y =Difference,label=Country), color = "black") + 
  # This will add USA labels to the USA points
  theme(legend.position = "none", 
  # This is useful for removing the legend
        axis.text.x = element_text(size = 15,angle = 30), 
        # this changes the size and angle of the x axis point labels 
        axis.text.y = element_text(size = 20), 
        axis.title.y = element_text(size =15), 
        axis.title.x = element_text(size =15), 
        strip.text.x = element_text(size = 15)) +
        # this changes the size of x axis labels for the facet
        ggtitle( "Change in BMI Over Time and Across Region Type and Gender")

        # Add a plot title  

We can now see that the rate of change from 1985 to 2017 appears to be larger in the women compared to the men in all regions. The group with the largest increase in the USA is the women living in rural areas.

Let’s check the difference with some statistical tests:

#first lets look at the normality of the difference in BMI
BMI_diff_long %>%
  ggplot(aes(x=Difference)) +
  geom_histogram() +
  facet_wrap(~ Sex + Type)

# interesting...the spread of difference is much broader for women than for men 
# some women in some countries especially need help with obesity
# while other countries are actually showing a loss of BMI
# these countries are the points on our previous plot with a difference below 0

#Let's see who:
BMI_diff_long %>% 
  arrange(Difference)

# A tibble: 1,200 x 4
   Country        Sex   Type                Difference
   <chr>          <chr> <chr>                    <dbl>
 1 Greece         Women Rural_difference        -1.5  
 2 Greece         Women National_difference     -1.3  
 3 Greece         Women Urban_difference        -1.20 
 4 Italy          Women Rural_difference        -1.10 
 5 Czech Republic Women Rural_difference        -0.8  
 6 Spain          Women Rural_difference        -0.800
 7 Lithuania      Women Urban_difference        -0.700
 8 Spain          Women Urban_difference        -0.700
 9 Lithuania      Women National_difference     -0.700
10 Spain          Women National_difference     -0.700
# … with 1,190 more rows

# In general, people appear to have lost weight in Greece

# Let's look at Q-Q plots:
BMI_diff_long %>%
  ggplot(aes(sample =Difference)) +
  stat_qq() + # to add the line we need to also include stat_qq_line()
  stat_qq_line() +
  facet_wrap(~ Sex + Type)

# ok so some of the data for Men looks fairly normal but not as much for Women

# Now we will check the variance

mood.test(pull(filter(BMI_diff_long, 
                      Sex == "Women", 
                      Type == "Rural_difference"), Difference), 
          pull(filter(BMI_diff_long, 
                      Sex == "Women", 
                      Type == "Urban_difference"), Difference))

    Mood two-sample test of scale

data:  pull(filter(BMI_diff_long, Sex == "Women", Type == "Rural_difference"),  and pull(filter(BMI_diff_long, Sex == "Women", Type == "Urban_difference"),     Difference) and     Difference)
Z = 3.5379, p-value = 0.0004033
alternative hypothesis: two.sided

#reject the null - conclude that variance is not equal

mood.test(pull(filter(BMI_diff_long, 
                      Sex == "Men", 
                      Type == "Rural_difference"), Difference), 
          pull(filter(BMI_diff_long, 
                      Sex == "Men", 
                      Type == "Urban_difference"), Difference))

    Mood two-sample test of scale

data:  pull(filter(BMI_diff_long, Sex == "Men", Type == "Rural_difference"),  and pull(filter(BMI_diff_long, Sex == "Men", Type == "Urban_difference"),     Difference) and     Difference)
Z = 0.59121, p-value = 0.5544
alternative hypothesis: two.sided

#accept the null - conclude that variance is equal

mood.test(pull(filter(BMI_diff_long, 
                      Sex == "Women", 
                      Type == "Rural_difference"), Difference), 
          pull(filter(BMI_diff_long, 
                      Sex == "Men", 
                      Type == "Rural_difference"), Difference))

    Mood two-sample test of scale

data:  pull(filter(BMI_diff_long, Sex == "Women", Type == "Rural_difference"),  and pull(filter(BMI_diff_long, Sex == "Men", Type == "Rural_difference"),     Difference) and     Difference)
Z = 10.843, p-value < 2.2e-16
alternative hypothesis: two.sided

#reject the null - conclude that variance is not equal


# let's compare the groups
# this is again still paired data
# this is because each value of "Rural_difference" and "Urban_difference" correspond to the same country

ks.test(pull(filter(BMI_diff_long, 
                    Sex == "Women",  
                    Type == "Rural_difference"), Difference), 
        pull(filter(BMI_diff_long,
                    Sex == "Women",  
                    Type == "Urban_difference"), Difference),
        paired = TRUE) 

    Two-sample Kolmogorov-Smirnov test

data:  pull(filter(BMI_diff_long, Sex == "Women", Type == "Rural_difference"),  and pull(filter(BMI_diff_long, Sex == "Women", Type == "Urban_difference"),     Difference) and     Difference)
D = 0.16534, p-value = 0.009043
alternative hypothesis: two-sided

ks.test(pull(filter(BMI_diff_long,
                    Sex == "Men",
                    Type == "Rural_difference"), Difference), 
        pull(filter(BMI_diff_long,
                    Sex == "Men", 
                    Type == "Urban_difference"), Difference),
        paired = TRUE)

    Two-sample Kolmogorov-Smirnov test

data:  pull(filter(BMI_diff_long, Sex == "Men", Type == "Rural_difference"),  and pull(filter(BMI_diff_long, Sex == "Men", Type == "Urban_difference"),     Difference) and     Difference)
D = 0.14586, p-value = 0.02995
alternative hypothesis: two-sided

ks.test(pull(filter(BMI_diff_long,
                    Sex == "Men", 
                    Type == "Rural_difference"), Difference), 
        pull(filter(BMI_diff_long, 
                    Sex == "Women",  
                    Type == "Rural_difference"), Difference),
        paired = TRUE)

    Two-sample Kolmogorov-Smirnov test

data:  pull(filter(BMI_diff_long, Sex == "Men", Type == "Rural_difference"),  and pull(filter(BMI_diff_long, Sex == "Women", Type == "Rural_difference"),     Difference) and     Difference)
D = 0.35204, p-value = 5.645e-11
alternative hypothesis: two-sided

Now we have performed seven comparisons (4 earlier) so we should apply our multiple testing correction

.05/7

[1] 0.007142857

We see that there is a significant difference in the change in BMI Rural communities between men and women. This change is larger for women.

Importantly- we noticed that it appears to be specific countries where BMI shows a particular increase especially for women. Which countries are those? how does that compare with the US? Clearly the US is among the countries with the highest differences.

BMI_diff_long %>% 
  filter(Country == "USA")

# A tibble: 6 x 4
  Country Sex   Type                Difference
  <chr>   <chr> <chr>                    <dbl>
1 USA     Men   Rural_difference          3.10
2 USA     Women Rural_difference          3.5 
3 USA     Men   Urban_difference          3.1 
4 USA     Women Urban_difference          3.40
5 USA     Men   National_difference       3.10
6 USA     Women National_difference       3.4 

In the US, focus should be placed on both urban and rural women to improve this public health issue.

Here we can see the countries have the largest differences in BMI from 1985-2017:

BMI_diff_long %>%
  filter(Difference>4.9) %>%
  arrange(-Difference)

# A tibble: 6 x 4
  Country     Sex   Type                Difference
  <chr>       <chr> <chr>                    <dbl>
1 Egypt       Women Rural_difference          5.9 
2 Honduras    Women National_difference       5.60
3 Honduras    Women Rural_difference          5.5 
4 Egypt       Women National_difference       5.10
5 Saint Lucia Women Rural_difference          5   
6 Honduras    Women Urban_difference          4.9 

#here we sorted by the highest to lowest values of BMI Difference

However it is important to see what the mean BMI values are were for these countries in 2017. It could be that the average was underweight in 1985… let’s take a look.

BMI_long %>% 
  filter(Country %in% pull(filter(BMI_diff_long, Difference>4.9), Country), 
                    Sex %in% pull(filter(BMI_diff_long, Difference>4.9), Sex),
                    Year =="2017") %>%
  arrange(-BMI)

# A tibble: 9 x 5
  Country     Sex   Region   Year    BMI
  <chr>       <chr> <chr>    <chr> <dbl>
1 Egypt       Women Urban    2017   32.3
2 Egypt       Women National 2017   31.7
3 Egypt       Women Rural    2017   31.3
4 Saint Lucia Women National 2017   30.5
5 Saint Lucia Women Rural    2017   30.5
6 Saint Lucia Women Urban    2017   30.5
7 Honduras    Women Urban    2017   28.3
8 Honduras    Women National 2017   27.7
9 Honduras    Women Rural    2017   26.9

#Here is the data just for the changes in BMI in Egypt
filter(BMI_diff_long, Country == "Egypt") %>% 
  arrange(-Difference)

# A tibble: 6 x 4
  Country Sex   Type                Difference
  <chr>   <chr> <chr>                    <dbl>
1 Egypt   Women Rural_difference          5.9 
2 Egypt   Women National_difference       5.10
3 Egypt   Women Urban_difference          4.20
4 Egypt   Men   Urban_difference          2.60
5 Egypt   Men   Rural_difference          2.5 
6 Egypt   Men   National_difference       2.5 

Thus rural women in Egypt had the greatest increase in BMI from 1985 to 2017 in this data (5.9) and the average BMI is now over the obesity threshold of 30.

The data suggests that rural women in Egypt and other countries may especially benefit from dietary resources and our interventions and programs to assist with weight management.

Now let’s make a plot that summarizes our findings. To do this we will simplify the other plots we made and then combine them together.

# simplified national means plot
Means_plot<-BMI_long %>% 
  filter(Sex %in% c("Men", "Women"), 
         Year %in% c("1985", "2017"), 
         Region == "National") %>%
ggplot(aes(x = Sex, y = BMI)) + 
  geom_boxplot(outlier.shape = NA, color = "black" , aes(fill = Sex)) + 
 scale_fill_manual(values=c("dodgerblue", "orchid2")) +
  facet_grid(~ Year) + 
  geom_hline(yintercept=30, linetype="dashed", color = "red", size =1) + 
  geom_jitter(data=BMI_long %>%
                filter(Sex %in% c("Men", "Women"), 
                       Year %in% c("1985", "2017"), 
                       Region == "National"),
                aes(x =Sex, y =BMI), 
                width = .2, size =2, shape =21, 
                color = "black", fill = "gray") +
  geom_jitter(data=subset(BMI_long %>%
                            filter(Sex %in% c("Men", "Women"), 
                                   Year %in% c("1985", "2017"), 
                                   Region == "National", 
                                   Country == "USA")), 
                             aes(x =Sex, y =BMI), 
                             width = .02, size =12, shape =21, 
                             color = "black", fill = "gray") + 
  # This will add points that are larger for the USA data
  geom_text(data=subset(BMI_long%>%
                          filter(Sex %in% c("Men", "Women"), 
                                 Year %in% c("1985", "2017"), 
                                 Region == "National", Country == "USA")),
                          aes(x =Sex, y =BMI,label=Country), 
                          color = "black") + 
  # This will add USA labels to the USA points
  theme_linedraw() +
  # This will make thebackground of the plot white and the facet labels black
  theme(legend.position = "none", 
  # This is useful for removing the legend
        axis.text.x = element_text(size = 15,angle = 30, vjust = 0.5), 
        # this changes the size and angle of the x axis point labels 
        axis.text.y = element_text(size = 20), 
        axis.title.y = element_text(size =15), 
        axis.title.x = element_text(size =15), 
        strip.text.x = element_text(size = 15))
   
#Simplified difference plot 
BMI_diff_long$Type <-BMI_diff_long$Type %>%
  str_replace( pattern = "_difference", replacement = "")

Diff_plot<-BMI_diff_long %>% 
  filter(Sex %in% c("Men", "Women"), 
         Type != "National") %>%
ggplot( aes(x = Type, y = Difference, col = Type)) + 
  geom_boxplot(outlier.shape = NA, color = "black" , aes(fill = Sex)) + 
  scale_fill_manual(values=c("dodgerblue", "orchid2")) +
  facet_grid(~ Sex) +
  geom_jitter(data=BMI_diff_long %>% 
    filter(Sex %in% c("Men", "Women"), 
           Type != "National"), 
    aes(x =Type, y =Difference), 
    width = .2, size =2, shape =21, 
    color = "black", fill = "gray") + 
  geom_jitter(data=subset(BMI_diff_long %>% 
    filter(Sex %in% c("Men", "Women"), 
           Type != "National", 
           Country == "USA")), 
    aes(x =Type, y =Difference), 
    width = .02, size =12, shape =21, 
    color = "black", fill = "gray") + 
  # This will add points that are larger for the USA data
  geom_text(data=subset(BMI_diff_long %>% 
    filter(Sex %in% c("Men", "Women"),
           Type != "National"), 
           Country == "USA"),
    aes(x =Type, y =Difference,label=Country), 
    color = "black") + 
  # This will add USA labels to the USA points
  theme_linedraw() +
    # This will make the background of the plot white and the facet labels black
  theme(legend.position = "none", 
  # This is useful for removing the legend
        axis.text.x = element_text(size = 15,angle = 30, vjust = 0.5), 
        # this changes the size and angle of the x axis point labels 
        axis.text.y = element_text(size = 20), 
        axis.title.y = element_text(size =15), 
        axis.title.x = element_text(size =15), 
        strip.text.x = element_text(size = 15))
  # this changes the size of x axis labels for the facet

#add labels
Diff_plot<-Diff_plot + 
        labs(title = "Change in BMI by region", 
             x = "", 
             y = "Change in BMI \n (1985 to 2017)") +
        theme(title = element_text (size = 12, face = "bold"))

Means_plot <-Means_plot + 
        labs(title = "Mean BMI over time", 
             x = "", 
             y = "Mean BMI") +
        theme(title = element_text (size = 12, face = "bold"))



obesity_text<-tibble(Year=c(1985),BMI=c(31),Sex=c("Men"),label=c("Obesity"))

Means_plot <-Means_plot + 
             geom_text(data = obesity_text,
                       label=pull(obesity_text,label), 
                       color = "red", 
                       aes( fontface ="bold.italic", size = 13))

cowplot::plot_grid(Means_plot, Diff_plot, labels = c("A", "B"))

Great, now we have put two plots together using the plot_grid() function of the cowplot package. This way we can clearly communicate two messages. The first being that BMI has increased over time globally and that many countries including the United States of America are approaching a mean BMI that is above the obesity threshold of 30. We can also see that women on average have larger BMI values than males, but that both genders show increased levels over time. In the second plot we can see that the increase in BMI is not just happening in urban communities, but in both rural and urban communities and particularly in women.

Summary

We have evaluated BMI average estimates from 200 different countries around the world. To do so we imported data from a pdf using the pdftools package. We used tidyverse packages such as dplyr, stringr, and tidy to clean the data and get it in a workable format. Our statistical analysis focused on evaluating differences in BMI in females around the world across time and between rural and urban areas. We found a significant difference both between years and between the type of community among women globally using t-tests and nonparametric tests. Thus BMI has increased in women since 1985. Although BMI estimates are significantly higher in Urban areas compared to rural areas, BMI estimates have increased in both regions. We learned that there are assumptions and considerations to keep in mind when performing tests that compare two groups.

We learned that the t-test relies on: 1) normality of the data for both groups (this is not as much of an issue if the number of observations is relatively large total n>30) 2) equal variance between the two groups (make sure you do the correct test if the data is not normal) 3) balanced sample sizes of the two groups 4) independent observations (or independent paired observations)

We learned that we can evaluate if our data is normally distributed by plotting the distribution and by creating Q-Q plots.

We learned that if our data is not normally distributed, we can consider these options: 1) We can still perform a t-test if our n is large 2) We can transform the data before performing a t-test 3) We can use a nonparametric test (Wilcoxon signed rank test, the Wilcoxon rank sum test, and the Two-sample Kolmogorov-Smirnov (KS) test) 4) We can perform a t-test with resampling methods (which should be especially considered when the groups are imbalanced)

We learned that if our groups do not have equal variance, that we need to use a modified t-test to account for this or the correct nonparametric test that does not rely on this assumption (Two-sample Kolmogorov-Smirnov test).

We also learned that if our groups are not well balanced, that we can consider resampling methods.

Finally its also important to remember that when we have paired data and intend to compare means, we should perform a paired test (paired t-test, Wilcoxon signed rank test, or paired KS test). Our data used in this case study was paired because observations were taken for the same countries for different categories of populations - thus, we wanted to compare these populations within the same country (male vs female, rural vs urban etc.). Other examples would be in a case-control study (if samples are matched for various demographics) or a study with repeated measures (for example, symptom measures from the same individual before and after a treatment).

Using the ggplot2 package we were able to visualize trends in the data. Importantly, the largest increase appears to be in the rural areas particularly for women. We were also able to identify how the US compared to other countries and which countries showed the largest increase in BMI over time. Analyses like this are important for defining which groups could benefit the most from interventions, education, and policy changes when attempting to mitigate public health challenges. You can see in the article that this data came from that many additional considerations would be involved to perform such an analysis.

Suggested Homework

Students can evaluate the change in BMI over time using the global data available for each year between 2015 and 2017. This data can be found here.